Econometric Reviews

ISSN: 0747-4938 (Print) 1532-4168 (Online) Journal homepage: [Link]/journals/lecr20

GMM inference in spatial autoregressive models

Süleyman Taşpınar, Osman Doğan & Wim P. M. Vijverberg

To cite this article: Süleyman Taşpınar, Osman Doğan & Wim P. M. Vijverberg (2018) GMM
inference in spatial autoregressive models, Econometric Reviews, 37:9, 931-954, DOI:
10.1080/00927872.2016.1178885

To link to this article: [Link]

View supplementary material

Published online: 11 Jul 2016.

Submit your article to this journal

Article views: 482

View related articles

View Crossmark data

Citing articles: 1 View citing articles

Full Terms & Conditions of access and use can be found at

[Link]
ECONOMETRIC REVIEWS
2018, VOL. 37, NO. 9, 931–954
[Link]

GMM inference in spatial autoregressive models

Süleyman Taşpınara , Osman Doğanb , and Wim P. M. Vijverbergc
a
Economics Program, Queens College, The City University of New York, New York, New York, USA; b Istanbul Ulasim A.S.,
Project Department, Istanbul, Turkey; c Ph.D. Program in Economics, The Graduate School and University Center, The City
University of New York, New York, New York, USA

ABSTRACT KEYWORDS
In this study, we investigate the finite sample properties of the optimal Asymptotic variance;
generalized method of moments estimator (OGMME) for a spatial econometric efficiency; GMM; inference;
model with a first-order spatial autoregressive process in the dependent vari- optimal GMME; SARAR;
spatial autoregressive
able and the disturbance term (for short SARAR(1, 1)). We show that the esti-
models; variance correction;
mated asymptotic standard errors for spatial autoregressive parameters can be Windmeijer correction
substantially smaller than their empirical counterparts. Hence, we extend
the finite sample variance correction methodology of Windmeijer (2005) JEL CLASSIFICATION
to the OGMME for the SARAR(1, 1) model. Results from simulation studies C13; C21; C31
indicate that the correction method improves the variance estimates in small
samples and leads to more accurate inference for the spatial autoregressive
parameters. For the same model, we compare the finite sample properties of
various test statistics for linear restrictions on autoregressive parameters. These
tests include the standard asymptotic Wald test based on various GMMEs,
a bootstrapped version of the Wald test, two versions of the C(α) test, the
standard Lagrange multiplier (LM) test, the minimum chi-square test (MC),
and two versions of the generalized method of moments (GMM) criterion test.
Finally, we study the finite sample properties of effects estimators that show
how changes in explanatory variables impact the dependent variable.

1. Introduction
The finite sample properties of generalized method of moments estimators (GMMEs) have been
extensively studied in the literature on nonspatial econometric models (in the 1996 special issue of
Journal of Business & Economic Statistics). For example, Hansen et al. (1996) show that the estimated
asymptotic variance of their continuous updating GMME for the asset pricing models are downward
biased in finite samples. For linear and nonlinear models of covariance structures, Altonji and Segal
(1996) and Clark (1996) show that the correlation between the optimal weight matrix and sampling
errors introduces finite sample bias and that the amount of bias is relatively greater when data have
long-tailed and skewed distributions. Furthermore, the optimal GMMEs (OGMMEs) of models of
covariance structures have smaller asymptotic standard errors, making inference and specification
testing problematic. Hall and Horowitz (1996); Brown and Newey (2002) confirm that the first-order
asymptotic theory often provides poor approximations to the distributions of test statistics obtained
from GMM estimators. Similarly, with regard to dynamic panel data models and panel count data
models, asymptotic standard errors for the OGMME may also be severely downward biased (more than
20%) relative to the empirical standard deviations (Arellano and Bond, 1991; Blundell and Bond, 1998;
Windmeijer, 2005; Bond and Windmeijer, 2005; Windmeijer, 2008).
GMM estimators for the spatial autoregressive models may suffer from similar deficiencies because
they are essentially multistep estimators. For instance, Kelejian and Prucha (1998, 2010) propose a

CONTACT Süleyman Taşpınar STaspinar@[Link] Economics Program, Queens College, The City University of New York,
New York, NY, USA.
Color versions of one or more of the figures in the article can be found online at [Link]/lecr.
© 2016 Taylor & Francis Group, LLC
932 S. TAŞPINAR ET AL.

multistep estimation method that involves a combination of IV and GMM estimation for the spatial
model with a spatial autoregressive process in the dependent variable and the disturbance term,
commonly referred to as SARAR(1, 1). In terms of estimation, the Kelejian and Prucha approach is
straightforward as computations do not require special techniques even for extremely large samples.
In principle, the GMM/IV estimator is inefficient relative to the maximum likelihood estimator (MLE)
(Prucha, 2014), but in finite samples the difference can be negligible (Das et al., 2003).
To increase the efficiency of GMMEs, Lee (2007a,b), Lin and Lee (2010), Liu et al. (2010), and Lee
and Liu (2010) suggest a set of linear and quadratic moment functions. The reduced forms of spatial
models specify expressions for spatial lags of the dependent variable (i.e., endogenous variables) in terms
of exogenous variables and the disturbance term. Then, linear moment functions may be constructed
on the basis of the deterministic part of the spatial lag terms, and quadratic moment functions may
be formulated from the stochastic part of the spatial lag variables. Both linear and quadratic quadratic
moment functions are chosen in such a way that the resulting GMME is asymptotically equivalent to
the MLE when disturbances are independent and identically distributed (i.i.d.) with a normal density.
When disturbances are simply i.i.d., Liu et al. (2010) and Lee and Liu (2010) show that the generalized
method of moments (GMM) estimator can be more efficient than the quasi MLE, respectively, for the
case of SARAR(1, 1) and SARAR(p, q).
The GMM estimation approach suggested in Lee (2007a,b), Lin and Lee (2010), and Liu et al. (2010)
requires initial consistent parameter estimates to obtain not only the optimal weight matrix but also the
linear and quadratic moment functions. The asymptotic theory in Lee (2007a,b), Lin and Lee (2010),
and Liu et al. (2010) assumes that these estimated matrices are infinitely precise and therefore does not
account for the variation stemming from the initial estimator used to construct these matrices. The
correction method suggested in Windmeijer (2005) specifically accounts for the presence of these initial
consistent parameter estimates in the optimal weight matrix. Therefore, it can provide improvements
for statistical inference. The aforementioned studies on the GMM estimation of spatial models indicate
that the GMMEs have desirable finite sample properties in terms of bias and root mean square error.
However, their performance in terms of estimated asymptotic standard errors and related test statistics
for statistical inference has not been explored. Hence, further research on the performance of the
estimated asymptotic standard errors of the GMME in spatial autoregressive models is warranted.
In this study, we evaluate the performance of the asymptotic approximations suggested for the GMM
estimators and MLE for the SARAR(1, 1) model. We extend the finite sample correction method of
Windmeijer (2005) to our model. The correction is based on a Taylor approximation of the first-order
conditions of the GMM objective function to account for the additional variation in the optimal weight
matrix. Therefore, the correction method accounts for the variation stemming from the initial GMM
estimation to obtain an estimate of the optimal weight matrix. In a Monte Carlo study, we compare
the estimated asymptotic standard errors and the corrected ones with the empirical standard deviations
under various scenarios. Our simulation results indicate that the estimates of the asymptotic standard
errors based on the OGMME are downward biased for the spatial autoregressive parameters. The amount
of bias is generally larger than 5%. For the correction method, the amount of bias is less than 5% in most
cases suggesting that the correction method can be useful for statistical inference.
For the SARAR(1, 1) specification, we consider various test statistics for testing linear restrictions on
spatial autoregressive parameters based on the GMM and ML methods. These tests include standard
asymptotic Wald tests based on the MLE and the GMME, a bootstrapped version of the Wald test, two
versions of the C(α) test, the standard LM test, the minimum chi-square test (MC), and two versions of
the GMM criterion test. We compare the performance of these tests through a Monte Carlo simulation
for the two null hypotheses that there is no spatial dependence (i) in the dependent variable or (ii) in the
disturbance term. Our results show that the standard asymptotic Wald test based on the OGMME tends
to over-reject, whereas the Wald test based on the corrected standard errors has proper size properties.
The GMM criterion tests and the C(α) test can be useful for testing the null hypothesis that there is no
spatial dependence in the disturbance terms. Our results also show that the C(α) test can be useful to
detect the presence of spatial dependence in the dependent variable.
 ECONOMETRIC REVIEWS 933

For the same specification, we study the finite sample properties of scalar measures suggested by
LeSage and Pace (2009) for the partial derivatives that show how changes in explanatory variables impact
the dependent variable. We consider (i) the simulation method of LeSage and Pace (2009) and (ii) the
Delta method to study the dispersions of the direct, indirect and total effects. Our results indicate that
both methods perform similarly and the dispersions calculated with the corrected standard errors for
the effects estimators can be useful as they have relatively smaller bias.
This article is organized in the following way. Section 2 presents the spatial autoregressive model
under consideration and discusses its assumptions. Section 3 reviews the GMM estimation approach and
lays out the details of the finite sample correction method for the SARAR(1, 1) specification. Section 4
reviews scalar measures suggested for the interpretation of parameters, and shows how their dispersions
can be calculated. Section 5 reviews the details of test statistics under consideration in the GMM
framework. Section 6 lays out the details of the Monte Carlo design and presents the results. Section 7
closes with concluding remarks. Some of the technical derivations and Monte Carlo results are relegated
to a web appendix.

2. The model specification and assumptions

The SARAR(1, 1) specification is given by
Yn = λ0 Wn Yn + Xn β0 + un , un = ρ0 Mn un + εn , (2.1)
where Yn is the n × 1 vector of a dependent variable, Xn is the n × k matrix of non-stochastic exogenous
variables with a conformable parameter vector β0 , Wn and Mn are the n × n spatial weight matrices of
known constants with zero diagonal elements, and εn is the n×1 vector of disturbances (or innovations).
Wn Yn and Mn un are known as, respectively, the spatial lag of the dependent variable and the spatial lag
of the disturbance term. λ0 and ρ0 are the spatial autoregressive parameters. Because spatial data is
characterized with triangular arrays, the variables in (2.1) have the subscript n (Kelejian and Prucha,
2010). Let 2 be the parameter space of the model. In order to distinguish the true parameter vector
from other possible values in 2, we state the model with the true parameter vector θ0 = (ρ0 , λ0 , β ′ 0 )′ .
Furthermore, for notational simplicity we let Sn (λ) = (In − λWn ), Rn (ρ) = (In − ρMn ), Gn (λ) =
Wn S−1 −1 S −1
n (λ), Hn (ρ) = Mn Rn (ρ), and Gn (ρ, λ) = Rn (ρ)Gn (λ)Rn (ρ). Also, at the true parameter values
(ρ0 , λ0 ), we denote Sn (λ0 ) = Sn , Rn (ρ0 ) = Rn , Gn (λ0 ) = Gn , Hn (ρ0 ) = Hn , and S
Gn = Rn Gn R−1
n .
The consistency and asymptotic normality of the GMME for (2.1) are established under the following
assumptions.1

Assumption 1. The elements εin of the disturbance term εn are distributed independently and identically
with mean zero and variance σ02 , and E |εin |4+ν < ∞ for some ν > 0 for all n and i.

Assumption 2. The spatial weight matrices Mn and Wn are uniformly bounded in row and column sums
in absolute value. Moreover, S−1 −1 −1 −1
n , Sn (λ), Rn , and Rn (ρ) exist and are uniformly bounded in row and
column sums in absolute value for all values of ρ and λ in a compact parameter space.

The regularity conditions in Assumptions 1 and 2 are motivated to restrict the spatial autocorrelation
in the model at a tractable level (Kelejian and Prucha, 1998). By this assumption, the third and fourth
moments of εin , denoted respectively by µ3 and µ4 , exist for all i and n. Assumption 2 also implies that
the model in (2.1) represents an equilibrium relation for the dependent variable. By this assumption,
the reduced form of the model can be written as Yn = S−1 −1 −1 2
n Xn β0 + Sn Rn εn . In the literature, weight
matrices are usually treated as exogenous and fixed. Lee (2004, 2007b) formulates the weight matrix as

1
For interpretations and implications of these assumptions, see Lee (2007a) and Kelejian and Prucha (2010).
2
The uniform boundedness of S−1 −1
n (λ) and Rn (ρ) is required for the ML estimator. In the case of GMM estimators, only the
−1 −1
uniform boundedness of Sn and Rn is required (Liu et al., 2010).
934 S. TAŞPINAR ET AL.

a function of the sample size, such that the sequence of weight matrices {Wn } is uniformly bounded
in both row and column sums and its elements wn,ij are O( h1n ). The sequence {hn } can be bounded or
divergent with the property that limn→0 hnn = 0, which implies that hn is allowed to diverge only at a
rate slower than that of n. For large group interactions for which limn→∞ hnn 6= 0, the estimators might
not be consistent. In this study, we assume that hn is bounded.
Throughout this study, the vector of moment functions
we consider for the GMM estimator takes
the form of g(θ0 ) = ε′ n P1n εn , . . . , ε′ n Pmn εn , ε′ n Qn ′ . Moment functions involving the n × n constant
matrices Pjn for j = 1, . . . , m are known as quadratic moment functions. The moment function Qn ′ εn is a
linear moment function where Qn is an n×r instrument matrix with r ≥ k+1 and has full column rank.
The matrices Pjn and Qn are chosen in such way that orthogonality conditions of population moment
functions are not violated. Let Pn be the class of n × n constant matrices with zero trace. The quadratic
moment functions involving matrices from Pn satisfy the orthogonality conditions when disturbance
terms are i.i.d. Assumption 4 states regularity conditions for these matrices. Assumption 5 characterizes
the parameter space.3

Assumption 3. The regressors matrix Xn is an n × k matrix consisting of uniformly bounded constant

elements. It has full column rank. Moreover, limn→∞ n1 X ′ n Xn exists and is nonsingular.

Assumption 4. Elements of the IV matrix Qn are uniformly bounded. Matrices Pjn for j = 1, . . . , m are
uniformly bounded in row and column sums in absolute value.

Assumption 5. The parameter space 2 is a compact subset of Rk+2 , and θ0 ∈ Int(2).

3. Spatial GMM and finite sample variance correction method

3.1. The moment conditions
Equation (2.1) implies Yn = S−1 −1
n (Xn β0 + Rn εn ). Thus, the endogenous variable Wn Yn on the right-
hand side of Eq. (2.1) may be written as Wn Yn = Wn S−1 −1
n (Xn β0 + Rn εn ) = Gn Xn β0 + Gn Rn εn .
−1

Equation (2.1) may also be written as Rn Yn = λ0 Rn Wn Yn + Rn Xn β0 + εn . This expression is the

basis for the choice of instruments (Lee, 2007a; Liu et al., 2010; Lee and Liu, 2010; Lin and Lee,
2010; Kelejian and Prucha, 2010): the linear moment matrix Qn is constructed from the expectation
of Zn = (Rn Wn Yn , Rn X
n ), which implies

that the orthogonality condition for the population moment
function holds: E Q′ n εn = Q′ n E εn = 0r×1 . Given consistent initial estimates of λ0 , ρ0 , and β0 , the
IV matrix Qn becomes available. In particular, Lee (2003) shows that the 2SLS estimator based on Qn is
the best estimator in the sense that its asymptotic covariance matrix is the smallest among the class of
2SLS estimators based on linear moment conditions.
Instruments that lead to quadratic moment functions follow from Assumption 1, in particular that
εin is i.i.d. with variance σ02 . For any Pjn ∈ Pn , Pjn εn is an
instrument′ that satisfies

the orthogonality


condition for a population moment function: E ε′ n Pjn ′ ε ′ ′ ′
n = E tr(ε n Pjn εn ) = tr Pjn E(εn ε n ) =


σ02 tr Pjn = 0 for j = 1, . . . , m, where tr(·) returns

the trace of an input

matrix. Moreover, when Pjn
is suitably chosen, we have E (Pjn εn )′ Rn Wn Yn = σ02 tr Pjn ′ R G R−1 6 = 0: P ε also instruments
n n n jn n
Rn Wn Yn .
The consistency of the GMM estimator does not depend on a particular Pjn , but the asymptotic
variance-covariance matrix is a function of utilized Pjn matrices. Therefore, the asymptotic efficiency of
the estimator is used as a criterion for the selection of Pjn matrices from Pn . Liu et al. (2010) suggest a set
of quadratic and linear moment matrices such that the resulting GMME is theoretically more efficient
than the quasi MLE. In this approach, the diagonal elements of the best quadratic and linear moment
matrices are also used to formulate additional moment functions to increase efficiency.

3
See Kelejian and Prucha (2010) for the specification of the parameter space of autoregressive parameters.
 ECONOMETRIC REVIEWS 935

To describe the best choice of matrices Pjn , we define the following notation. For the regressors matrix,
S
let Xn = Rn Xn = (Rn ln , Rn Xn⋆ ), where ln denotes n × 1 vector of ones and Xn⋆ is the submatrix of Xn with
the intercept column deleted. The submatrix Rn Xn⋆ = S Xn⋆ is assumed to be n × k⋆ matrix with column
S⋆ ⋆ 4
vectors denoted by Xnl for l = 1, . . . k . Let D(·) be the operator which either creates a matrix from the
diagonal elements of an input matrix or returns a diagonal matrix if the input is a vector. Furthermore, let
vec(·) be the operator that creates a column vector from the elements of an input matrix, and let vecD (·)
be the operator that creates a column vector from the diagonal elements of an input matrix. Finally, for
(t) (s)
any n × n matrix An , let An = An − n1 tr(An )In , and let An = An + A′ n .
⋆ =S (t) ⋆ (t)
Define the following quadratic moment matrices: P1n Gn , P2n = D(S Gn ) , P3n ⋆ = D(S
Gn S
Xn β0 )(t) ,
(t) (t)
P4n = Hn , P5n = D(Hn ), and Pl+5,n = D(S
⋆ ⋆ ⋆ ⋆ (t) ⋆
Xnl ) for l = 1, . . . , k . Note that Pjn ∈ Pn ⋆


for j = 1, . . . , k + 5. The linear moment matrices are Q⋆n = Q⋆1n , Q⋆2n , Q⋆3n , Q⋆4n , Q⋆5n =
⋆
(t)

S
Xn⋆ , S
Gn S
Xn β0 , ln , vecD (S G(t)
n ), vecD (Hn ) , which, one might note, expands on the matrix Qn given
above. Then, as shown by Liu et al. (2010), the best set of moment functions is given by


gn⋆ (θ ) = ε′ n (θ )P1n
⋆
εn (θ ), . . . , ε′ n (θ )Pk⋆⋆ +5,n εn (θ ), ε′ n (θ )Q⋆n ′ ,
where εn (θ ) = Rn (ρ)Sn (λ)Yn − Rn (ρ)Xn β. This set of moment functions is best in the sense that any
other moment function that can be added to this set does not increase the asymptotic efficiency of the
GMME. h i h ⋆ i
′ ∂gn (θ0 )
Let ⋆n = E gn⋆ (θ0 )gn⋆ (θ0 ) and 8⋆n = −E ∂θ ′ . It can be shown that5
   
⋆ (µ4 − 3σ04 )ϑnm ′ ϑnm µ3 ϑnm ′ Q⋆n 4 4nm 0m×r
n = ′ + σ0 ′ , (3.1)
µ3 Q⋆n ϑnm 0r×r 0r×m σ0−2 Qn⋆ Q⋆n
 ⋆(s) ⋆(s)  
where m = k⋆ + 5, r = k⋆ + 4, 4nm = vec(P1n ), . . . , vec(Pk⋆ +5,n ) ′ vec(P1n ⋆ ), . . . , vec(P⋆
k⋆ +5,n ) , and
 ⋆ ), . . . , vec (P⋆

ϑnm = vecD (P1n D k⋆ +5,n ) . The negative expectation of the gradient matrix is given by
 ⋆(s)
⋆(s)


σ02 tr H ′ n P1n σ02 tr SG′ n P1n 01×k
 .. .. 
 . . 
⋆
8n = −  

. (3.2)
⋆(s) ⋆(s) 
 σ02 tr H ′ n Pk⋆ +5,n σ02 tr S G′ n Pmn 01×k 
′ ′
0r×1 Q⋆n S Gn SXn β0 Q⋆n S Xn
The results in (3.1) and (3.2) indicate that consistent estimates of ⋆n and 8⋆n require consistent estimates
of θ0 , σ02 , µ3 , and µ4 . In practice, an initial estimator of θ0 can be used to recover consistent estimates of
these parameters.6 Let θ̂1 be an initial GMME (IGMME). Let n (θ̂1 ) denote the estimate of ⋆n recovered
from θ̂1 . Then, the OGMME is defined by
′
θ̂n⋆ = argminθ∈2 gn⋆ (θ ) −1 ⋆
n (θ̂1 ) gn (θ ). (3.3)
Under Assumptions 1–5, Liu et al. (2010) show that the GMM estimator defined in (3.3) is consistent
and asymptotically normally distributed, namely7
  !
√ ⋆
d 1 ⋆′ ⋆−1 ⋆ −1
n θ̂n − θ0 −→ N 0(k+2)×1 , lim 8n n 8n . (3.4)
n→∞ n

4
In the case where Xn has no column of ones, Rn Xn = Rn Xn⋆ = S Xn⋆ and k⋆ = k.
5
Lemma 1 in Web Appendix A can be used to derive ⋆n matrices in this section.
6
An initial estimator of θ0 can be the one suggested in Prucha (2014). Another candidate is the initial GMM estimator (IGMME)
suggested at the end of this section. In addition, the moment matrices also are functions of unknown parameters. With initial
estimates, these matrices also became available.
7
The identification of parameters in the GMM framework requires limn→∞ n1 E(gn (θ0 )) = 0 (Newey and McFadden, 1994).
Lee (2007a) investigates this condition and state the identification conditions for parameters. Here, we simply assume that
the parameters are identified.
936 S. TAŞPINAR ET AL.

√

An estimate for the covariance of n θ̂n⋆ − θ0 is needed for statistical inference. The result in (3.4) sug-
 ′ 
⋆ ⋆ −1 , where 8⋆ (θ̂ ⋆ ) is an estimate
gests that an estimate may be obtained from n1 8⋆n (θ̂n⋆ ) −1 n (θ̂1 ) 8n (θ̂n ) n n
 ′ 
⋆ ⋆ −1 ,
of 8⋆n recovered from θ̂n⋆ . Alternatively, one might form an estimate from n1 8⋆n (θ̂n⋆ ) −1 ⋆
n (θ̂n ) 8n (θ̂n )
where n (θ̂n⋆ ) is an estimate of ⋆n recovered from θ̂n⋆ . Although these two procedures are asymptotically
equivalent, little seems to be known about the relative merits of these two procedures in small samples
(Newey and McFadden, 1994).
The estimator defined in (3.3) is the best GMM estimator in the class of OGMME formulated from
the set of linear and quadratic moment functions. When disturbances are i.i.d. normal, the best set of
(t) (t)

moment functions is given by P1n = S Gn , P2n = Hn , and Qn = S Gn SXn β0 , S Xn . Liu et al. (2010) show
that any other moment function that can be added to this set is redundant. Furthermore, the MLE is
also characterized by this set of moment functions, which implies that the GMME based on this set
is asymptotically equivalent to the MLE. When disturbances are simply i.i.d., the additional moment
(t)

(t) (t)
(t)
(t)


⋆ (t) for
functions formulated from D S Gn , D S Gn S
Xn β0 , D Hn , vecD S Gn , vecD Hn , and D S Xnl
l = 1, . . . , k⋆ increase asymptotic efficiency. The resulting GMME based on the bigger set gn⋆ (θ ) is
asymptotically more efficient than the quasi MLE (Liu et al., 2010). Note that, the efficiency gain from
these additional moment functions depends on the specification of weight matrices and regressors.
For the case of SARAR(1,
0), the best set of moment functions follows from the more general case
above. Let δ0 = λ0 , β ′ 0 ′ be the true parameter vector. Then, the best set of moments gn⋆ (δ) uses
⋆ (t) ⋆ (t) ⋆ ⋆ ⋆ )(t) for l = 1, . . . , k⋆ , and
P1n = Gn , P2n = D(Gn ), P3n = D(Gn Xn β0 )(t) , Pj+3,n = D(Xnl
(t)

Q⋆n = (Q⋆1n , Q⋆2n , Q⋆3n , Q⋆4n ) = Xn⋆ , Gn Xn β0 , ln , vecD (Gn ) , where εn (δ) = Sn (λ)Yn − Xn β. Similarly,
the best set of the moment functions for the case of SARAR(0, 1) may be obtained. Let η0 = (ρ0 , β ′ 0 )′
⋆ = H (t) , P⋆ = D H (t) ,


be the true parameter vector. Then, the best set of moments gn⋆ (η) uses P1n n 2n n


⋆ (t) for l = 1, . . . , k⋆ , and Q⋆ = Q⋆ , Q⋆ , Q⋆

(t)

⋆
Pl+2,n = DS Xnl n 1n 2n 3n = S Xn⋆ , ln , vecD (Hn ) where
εn (η) = Rn (ρ)Yn − Rn (ρ)Xn β.
Note that Pjn⋆ and Q⋆ also involve unknown parameters. In practice, an initial consistent estimator is
n
used to construct these matrices. For example, an IGMME of θ0 based on quadratic moment matrices



(t)
(t)
P1n = Wn ′ Wn − D Wn ′ Wn , P2n = Mn ′ Mn − D Mn ′ Mn , P3n = W ′

n Wn , and P4n = Mn Mn ,
′

and the linear moment matrix Qn = Wn Mn Xn , Wn Xn , Mn Xn , Xn may be used to construct these

matrices when disturbance terms satisfy Assumption 1. Hence, the optimal GMM estimator defined
in (3.4) becomes feasible when there is an initial consistent estimator of θ0 . The optimal weight, the
quadratic, and linear moment matrices are constructed from the initial estimator. The replacement of
⋆ , Q⋆ , and ⋆ with their consistent estimates does not affect the asymptotic properties of the optimal
Pjn n n
GMM estimators.8 However, it may affect the small-sample behavior of the estimators.

3.2. Finite sample variance correction method

In this section, the finite variance correction methodology of Windmeijer (2005) is extended for our
spatial autoregressive model. We provide a general argument by
considering the general set of population
moment functions gn (θ0 ) = ε′ n P1n εn , . . . , ε′ n Pmn εn , ε′ n Qn ′ , where Qn is an n × r matrix of linear
instruments, and Pjn ∈ Pn for j = 1, . . . , m.
The extension of the finite sample variance correction method to our case is not straightforward for
mainly two reasons. First, our generic set of moment functions cannot be written at the observational
level in such a way that the resulting moment functions are independent across observations. This arises
simply because of the spatial dependence. For the framework considered in Windmeijer (2005), the
set of moment functions is defined at the observational level such that gin (θ ) is the sample moment
function for the ith observation. Within this framework, a good candidate for an estimate of n is

8
Lee (2007a, Proposition 2) and Liu et al. (2010, Proposition 2, p. 306) give the asymptotic arguments.
 ECONOMETRIC REVIEWS 937

P
given by n1 ni=1 gin (θ̂n )g ′ in (θ̂n ), where θ̂n is a consistent initial estimator (Newey and McFadden, 1994,
p. 2155). This candidate is no longer valid for our case because of the spatial dependence implied by
the model specification. The second reason is a byproduct of the first and is related to the feasibility
of the variance
 correction
 formula for our spatial model. In our case, an estimate of n is recovered
from E gn (θ0 )g ′ n (θ0 ) . The correction method is based on the idea of taking a first-order Taylor
approximation of n (θ̂1 ) around the true parameter vector. As we will show below, this step is not as
straightforward as the one stated in Windmeijer (2005). In what follows, we present the details on how
the correction method is extended to our model.
Let 9n−1 be an arbitrary non-stochastic weighting matrix for the GMM objective function. The
weighting matrix plays the role of a metric by which the sample moment functions are made as close
as possible to zero. Furthermore, assume that 9n converges to a constant positive definite matrix
9. The IGMM estimator is obtained by minimizing the objective function with respect to θ . Let
Q9n (θ ) = g ′ n (θ )9n−1 gn (θ ) be the objective function when the weight matrix is 9n . Then, the initial
√ d
estimator is defined

as θ̂1 = argminθ Q9n (θ ).9 Under our stated assumptions, we have n(θ̂1 − θ0 ) − →
N 0(k+2)×1 , ϒ , where
 −1       −1
1 ′ −1 1 1 ′ −1 1 −1 1 1 ′ −1 1
ϒ = lim 8n 9n 8n 8n 9n n 9n 8n 8n 9n 8n . (3.5)
n→∞ n n n n n n n
Once we have a consistent estimator, we can construct an estimate of n for the optimal GMM
estimation. Let n (θ̂1 ) be the estimate of n constructed from θ̂1 . Similarly, we will denote the unknown
n with n (θ0 ) in the rest of this section. Then the OGMME is defined as θ̂2 = argminθ Qn (θ̂1 ) .
∂g (θ)
Let Cn (θ ) = ∂θ n
′ be the matrix of first derivatives and Fn (θ ) = ∂C∂θ
n (θ)
be the matrix of second
derivatives.10 The first and second derivatives of the GMM objective function are given by

1 ∂Q9n
B9n (θ0 ) = = Cn ′ (θ0 )9n−1 gn (θ0 ), (3.6)
2 ∂θ θ0

1 ∂ 2 Q9n  
A9n (θ0 ) = = Cn ′ (θ0 )9n−1 Cn (θ0 ) + Fn ′ (θ0 ) Ik+2 ⊗ 9n−1 gn (θ0 ) , (3.7)
2 ∂θ ∂θ ′ θ0

where B9n (θ0 ) is a (k + 2) × 1 vector and A9n (θ0 ) is a (k + 2) × (k + 2) matrix. When the weight matrix
is n in the objective function, these matrices are denoted by Bn (θ0 ) and An (θ0 ), respectively. From
the definition of the GMM objective function, we have Bn (θ̂1 ) (θ̂2 ) = 0(k+2)×1 and B9n (θ̂1 ) = 0(k+2)×1 .
Then, ignoring the remainder terms, a first-order Taylor series approximation of (3.6) around θ0 for the
OGMME and IGMME yields, respectively,


θ̂2 − θ0 = −A−1 (θ0 )Bn (θ̂1 ) (θ0 ), (3.8)
n (θ̂1 )


θ̂1 − θ0 = −A−1
9n (θ0 )B9n (θ0 ). (3.9)

Note that the expansion in (3.8) is also a function of θ̂1 . Therefore, a further expansion of θ̂1 around θ0
yields



θ̂2 − θ0 = −A−1 n (θ0 ) (θ0 )Bn (θ0 ) (θ0 ) + Tn (θ0 ) (θ0 ) θ̂1 − θ0 , (3.10)

9
For notional simplicity, we suppress the subscript n for the estimators. In order to make a difference between the IGMME
and the OGMME, we use θ̂1 for the IGMME and θ̂2 for the OGMME in this section.
10
In Web Appendix B, we provide the explicit expressions for Cn (θ ) and Fn (θ ) for our generic set of moment functions.
938 S. TAŞPINAR ET AL.

 
∂ A−1
n (θ) (θ0 )Bn (θ) (θ0 )
where Tn (θ0 ) (θ0 ) = ∂θ ′ is a (k+2)×(k+2) matrix. The jth column of Tn (θ0 ) (θ0 )
θ0
is given by (see Windmeijer, 2005, p. 28)
(
  −1 ∂n (θ )
Tn (θ0 ) (θ0 ) •, j = −An (θ0 ) (θ0 ) C′ n (θ0 )−1
n (θ0 ) −1
n (θ0 )Cn (θ0 )
∂θj θ0
!)
∂n (θ )
+F ′ n (θ0 ) Ik+2 ⊗ −1
n (θ0 ) −1
n (θ0 )gn (θ0 ) A−1
n (θ0 ) (θ0 )Bn (θ0 ) (θ0 )
∂θj θ0
∂n (θ )
+ A−1 ′ −1
n (θ0 ) (θ0 )Cn (θ0 )n (θ0 ) −1
n (θ0 )gn (θ0 ). (3.11)
∂θj θ0

Notice that the first term in Tn (θ0 ) (θ0 )[•, j] is a function of A−1
n (θ0 ) (θ0 )Bn (θ0 ) (θ0 ). This result indicates
−1
that the term An (θ0 ) (θ0 )Bn (θ0 ) (θ0 ) is the bias of an infeasible GMM estimator based on the unknown
 
optimal weight matrix n (θ0 ) = E gn (θ0 )g ′ n (θ0 ) . Hence, the first term is a function of this bias, which
tends to be small and generally remains constant as the number of moment functions increases. The
remaining term in Tn (θ0 ) (θ0 )[•, j] is the dominant term, and it becomes larger as the number of moment
functions increases (Windmeijer, 2005, p. 28).

The result in (3.10) may be used to formulate a variance formula for θ̂2 . The term Tn (θ0 ) (θ0 ) θ̂1 −θ0
in (3.10) has an order of Op ( n1 ) and vanishes as n → ∞, but it does not equal zero in finite samples.
Therefore, accounting for this term in a variance formula for θ̂2 may improve inference in finite samples.
Equations (3.9) and (3.10) result in the following finite sample variance θ̂2 :
Varc (θ̂2 ) = A−1 ′ −1 −1
n (θ0 ) (θ0 )C n (θ0 )n (θ0 )Cn (θ0 )An (θ0 ) (θ0 )
+ Tn (θ0 ) (θ0 )A−1 ′ −1 −1
9n (θ0 )C n (θ0 )9n Cn (θ0 )An (θ0 ) (θ0 )
+ A−1 ′ −1 −1 ′
n (θ0 ) (θ0 )Cn (θ0 )9n Cn (θ0 )A9n (θ0 )T n (θ0 ) (θ0 )
+ Tn (θ0 ) (θ0 )Var(θ̂1 )T ′ n (θ0 ) (θ0 ). (3.12)
The variance formula in (3.12) is infeasible as it involves the unknown θ0 . Windmeijer (2005) circum-
vents this problem by simply replacing the unknown parameters with their consistent estimates. The
same approach cannot be directly adopted for our case as the functional form for the estimate of n is
now different. An estimate for the unknown term Tn (θ0 ) (θ0 )[•, j] is recovered from the estimates of θ̂1
and θ̂2 . The first term in Tn (θ0 ) (θ0 )[•, j] vanishes as Bn (θ0 ) (θ̂n ) = 0(k+2)×1 for any consistent estimator
θ̂n . Hence, the first term of (3.11) does not play any role in the calculation of the variance of θ̂2 . Then, a
consistent estimate of the jth column of Tn (θ0 ) (θ0 ) for j = 1, . . . , k + 2 is given by
  ∂n (θ )
Tn (θ̂1 ) (θ̂2 ) •, j = A−1 (θ̂2 )Cn ′ (θ̂2 )−1 ˆ
n (θ1 ) −1
n (θ̂1 )gn (θ̂2 ). (3.13)
n (θ̂1 ) ∂θj θ̂2
More explicitly, the terms in (3.13) are given by
 
1 ∂n (θ ) 1 ∂gn (θ ) ′ ∂g ′ n (θ )
= E g n (θ ) θ0
+ gn (θ ) + op (1), (3.14)
n ∂θj θ̂2 n ∂θj ∂θj θ0
h i
An (θ̂1 ) (θ̂2 ) = C′ n (θ̂2 )−1 ′ −1
n (θ̂1 )Cn (θ̂2 ) + Fn (θ̂2 ) Ik+2 ⊗ n (θ̂1 )gn (θ̂2 ) , (3.15)
∂gn (θ) ∂Cn (θ)
where Cn (θ̂2 ) = ∂θ ′ θ̂2 and Fn (θ̂2 ) = ∂θ θ̂2
. Note that the estimates Cn (θ̂2 ), Fn (θ̂2 ), and gn (θ̂2 )
are simply obtained by evaluating Cn (θ ), Fn (θ ), and gn (θ ) at the consistent estimate θ̂2 . On the other
 ECONOMETRIC REVIEWS 939

hand, the consistent estimates n (θ̂1 ) and ∂∂θn (θ)

j θ̂2
are recovered from the corresponding results derived
under the expectation operator.
 ∂gn (θ) ′ ∂g ′ (θ) 
Next, we provide explicit expressions for E ∂θ j
g n (θ ) θ + gn (θ ) ∂θn j θ for j = 1, . . . , k + 2.
0 0
Let bk(·) be the matrix operator that creates a generalized block diagonal matrix from a given list of
 
 Ar×k 0r×l 
matrices. For example, if c = Ar×k , Br×l , then bk c = 0r×k Br×l . With this new notation, we have
the following result with respect to λ:
 
∂gn (θ ) ′ ∂g ′ n (θ )
E g n (θ )
θ0
+ gn (θ )
∂λ ∂λ θ0
"  (s) 

 #
(µ − 3σ 4) c ′c µ bk c ′ l ⊗Q +c ′S
G ′Q
=−  4 0 1nm 2nm

 3 1nm m n 2nm n n
µ3 Qn ′S
Gn c2nm + lm ′ ⊗ Qn ′ bk c1nm 0r×r
   
−2
′


[c ](s)
σ I ⊗ S
G S
X β bk d′ 1nm lm ⊗ Qn
 3nm 0 m n n 0 
− σ04 
 −2 ′





 

(s)
′
σ0 lm ⊗ Qn bk d1nm Im ⊗ Gn Xn β0 S S σ0−2 Qn ′S Gn Qn
   (s) 

′ ′
S S
−  µ3 Im ⊗ Gn Xn β0 d 1nm c2nm 0m×r  , (3.16)
0r×m 0r×r
 (s) (s)   
where c1nm = vecD (P1n S Gn ), . . . , vecD (PmnS Gn ) , c2nm = vecD (P1n ), . . . , vecD (Pmn ) , c3nm =
 (s) (s)  (s) (s)   (s) (s) 
vec(P1n SGn ), . . . , vec(PmnS Gn ) ′ vec(P1n ), . . . , vec(Pmn ) , d1nm = P1n , . . . , Pmn , Im is the m × m
identity matrix, lm is the m×1 vector of ones, and ⊗ denotes the kronecker product. In a similar fashion,
with respect to ρ we obtain
 
∂gn (θ ) ′ ∂g ′ n (θ )
E g n (θ ) θ + gn (θ )
∂ρ 0 ∂ρ θ0
  (s)    " #
4) d ′c ′+c ′H ′ Q (s)
(µ 4 − 3σ 0 2nm 2nm µ 3 d2nm 2nm n n 4 [d3nm ] 0m×r
=− − σ0 (s) ,
µ3 Qn ′ [Hn c2nm + d2nm ] 0r×r 0r×m σ0−2 Qn ′ Hn Qn
(3.17)
 (s) (s)   (s) (s)  (s)
where d2nm = vecD (P1n Hn ), . . . , vecD (Pmn Hn ) , and d3nm = vec(P1n Hn ), . . . , vec(Pmn Hn ) ′ vec(P1n ),
(s) 
. . . , vec(Pmn ) . Finally, with respect to βj for j = 1, . . . , k, we obtain
 
∂gn (θ ) ′ ∂g ′ n (θ )
E g n (θ ) θ + gn (θ )
∂βj 0 ∂βj θ0
   (s)   

′ ′ 2

′ ′



 µ3 Im ⊗ Rn Xnj d 1nm c2nm σ0 Im ⊗ Rn Xnj bk d 1nm lm ⊗ Qn 
= −    ,

2 ′ ′




σ0 lm ⊗ Qn bk d1nm Im ⊗ Rn Xnj 0r×r

(3.18)
where Xnj is the jth column of Xn . Note that under our stated assumptions, estimates of (3.16), (3.17), and
(3.18) are recovered by replacing (θ ′ 0 , σ02 , µ3 , µ4 )′ with their consistent estimates (θ̂ ′ n , σ̂n2 , µ̂3 , µ̂4 )′ .11

11
The asymptotic argument is similar to the one given in the proof of Proposition 2 of Liu et al. (2010). Hence, it is omitted.
940 S. TAŞPINAR ET AL.

 ∂gn (θ) ′
When θ̂2 is used to estimate (3.16), (3.17), and (3.18), the resulting estimate of E ∂θ j
g n (θ ) θ +
0
∂g ′ n (θ)  ∂n (θ)
gn (θ ) ∂θj θ is denoted by ∂θj θ̂ . Then, inserting these results into (3.12) yields an estimate of
0 2
the variance of θ̂2 that incorporates the variation in the IGMME θ̂1 :
c c (θ̂2 ) = A−1
Var (θ̂2 )C′ n (θ̂2 )−1
n (θ̂1 )Cn (θ̂2 )A
−1
(θ̂2 )
n (θ̂1 ) n (θ̂1 )

+ Tn (θ̂1 ) (θ̂2 )A−1 ′ −1

9n (θ̂1 )C n (θ̂1 )9n Cn (θ̂2 )A
−1
(θ̂2 )
n (θ̂1 )

+ A−1 (θ̂2 )Cn ′ (θ̂2 )9n−1 Cn (θ̂1 )A−1 ′

9n (θ̂1 )T n (θ̂1 ) (θ̂2 )
n (θ̂1 )
c θ̂1 )T ′
+ Tn (θ̂1 ) (θ̂2 )Var( (3.19)
n (θ̂1 ) (θ̂2 ).

Using (3.9), the estimate of variance of θ̂1 is obtained as

c θ̂1 ) = A−1 (θ̂1 )Cn ′ (θ̂1 )9n−1 n (θ̂1 )9n−1 Cn (θ̂1 )A−1 (θ̂1 ).
Var( (3.20)
9n 9n

The estimate of the variance of θ̂2 in (3.19) accounts for the variation that stems from the initial
estimator θ̂1 . Windmeijer (2005, 2008) shows that this finite sample correction provides a better
approximation to the empirical standard deviations of the estimator when the estimator itself is not
substantially biased. The Monte Carlo studies in Lee (2007a), Lin and Lee (2010), and Liu et al. (2010)
indicate that the bias of the best GMME of the spatial autoregressive model is negligible. Hence, the
correction formula provided in (3.19) can be useful for statistical inference.

4. Direct, indirect and total effects

In this section, we describe the marginal effects (impact measures) accounting within the context of our
spatial model. The partial derivative of the dependent variable with respect to kth explanatory variable
yields the following matrix of marginal effects:

∂Yn
= S−1
n βk0 , (4.1)
∂X ′ nk
where βk0 is the kth component of β0 . The result in (4.1) indicates that the marginal effects of
a change
in Xnk is represented by an n × n matrix. The diagonal elements of this matrix ∂Yni /∂Xn,ki contain
the own-partial

derivatives, while the off-diagonal elements represent the cross-partial derivatives
∂Ynj /∂Xn,ki . LeSage and Pace (2009) define the average of the main diagonal elements of this matrix as
a scalar summary measure of direct effects, and the average of off-diagonal elements as a scalar summary
measure of indirect effects. The sum of direct and indirect effects is labeled as the total effects.
We consider two methods for the calculation of dispersions of these impact measures: (i) a simulation
approach suggested by LeSage and Pace (2009), and (ii) the Delta method (Debarsy et al., 2015). The
simulation approach utilizes the parameter estimates and the estimated covariance matrix of a consistent

estimator. Let Ln be a lower-triangular matrix recovered from the Cholesky decomposition of Var θ̂n ,
and ϑn be a random vector that has a multivariate standard normal distribution. Then, random draws
of the parameter vector are generated according to

θ r = Ln ϑnr + θ̂n , for r = 1, . . . , R. (4.2)

For each θ r , the direct, indirect, and total effects can be calculated using (4.1). The mean and the standard
deviation calculated from each sequence of impact measures can be used as the point estimate and the
standard error of the corresponding impact measure.
 ECONOMETRIC REVIEWS 941

The second approach is the conventional method based
on the Delta method. The result in (4.1)
indicates that the estimate of direct effect is n1 tr S−1
n (λ̂n )β̂nk . By the mean value theorem,
1 h  
i √  
√ tr S−1 −1
n (λ̂n )β̂nk − tr Sn βk0 = A1n × n λ̂n − λ0 , β̂nk − βk0 ′ + op (1)
n
d
 
→ N 0, lim A1n Bn A′ 1n ,
− (4.3)
n→∞
1 1



−1 , and B is the asymptotic covariance of √n λ̂ − λ , β̂ −
where A1n = n tr S−1 n Gn βk0 , n tr Sn n n 0 nk

′
βk0 . The result in (4.3) indicates that the asymptotic variance of direct effects can be estimated

1

by n1 Â1n B̂n Â′ 1n , where Â1n = n1 tr S−1 −1
n (λ̂n )G
n (λ̂n )β̂nk , n tr Sn (λ̂n ) , and B̂n is the estimated
√
asymptotic covariance of n λ̂n − λ0 , β̂nk − βk0 ′ .
Applying the mean value theorem to the estimate of total effects n1 β̂nk l′ n S−1 n (λ̂n )ln yields
1 h i √  
√ β̂nk l′ n S−1 ′ −1
n (λ̂n )ln − βk0 l n Sn ln = A2n × n λ̂n − λ0 , β̂nk − βk0
′
+ op (1)
n
d
 
→ N 0, lim A2n Bn A′ 2n ,
− (4.4)
n→∞
1 1 ′ −1
 1


where A2n = n βk0 n S−1 l′ ′ −1
n Gn ln , n l n Sn ln . Hence, Var n β̂nk ln Sn (λ̂n )ln can be estimated by
1 ′
 1 ′ −1 1 ′ −1
n Â2n B̂n  2n , where Â2n = n β̂nk l n Sn (λ̂n )Gn (λ̂n )ln , n l n Sn (λ̂)ln .

1
The estimate of indirect effects is given by n β̂nk l n Sn (λ̂n )ln − tr S−1
′ −1
n (λ̂n )β̂nk . The results in (4.3)
and (4.4) imply that
1 h  

i
√ β̂nk l′ n S−1 −1
n (λ̂n )ln − tr Sn (λ̂n )β̂nk − βk0 l′ n S−1 −1
n ln − tr Sn )βk0
n
√   d
 
= (A2n − A1n ) × n λ̂n − λ0 , β̂nk − βk0 ′ + op (1) − → N 0, lim (A2n − A1n ) Bn (A2n − A1n ) ′ .
n→∞
(4.5)
1


1



Hence, an estimate of Var n β̂nk l n Sn (λ̂n )ln −tr Sn (λ̂n )β̂nk is given by n Â2n −Â1n B̂n Â2n −Â1n ′ .
′ −1 −1

5. Hypothesis testing with the GMME

In this section, we consider a Wald test procedure that is familiar within the general GMM framework
for testing linear and nonlinear restrictions on the parameter vector θ0 . We provide details about how to
adjust the formulation of this test statistic for the spatial SARAR(1, 1) model.
Let r : Rp → Rq be a twice continuously differentiable function, and assume that its derivative matrix
R = ∂r(θ )/∂θ ′ has rank q.12 Let θ̂n be an unrestricted optimal GMME. Then, the Wald test statistic for q
linear restrictions of the form r(θ0 ) = 0 has a limiting chi-square distribution with q degrees of freedom,
namely
 −1 d
Wald = r̂n ′ R̂′ n Var(θ̂n )R̂n → χq2 ,
r̂n − (5.1)

where r̂n = r(θ̂n ) and R̂n = ∂r(θ)

∂θ ′ θ̂n . The Wald test statistic may be calculated with the variance formulas
given (3.4), (3.5), and (3.19).
We also consider a bootstrapping method for the distribution of the Wald test statistic in our GMM
framework for the null hypotheses of H0 : λ0 = 0 and H0 : ρ0 = 0. Hall and Horowitz (1996), Brown and
Newey (2002), and Bond and Windmeijer (2005) outline bootstrap methods to determine finite sampling
distribution of Wald statistics. In the spatial econometrics literature, examples of these methods may be

12
Note that for our spatial model, p = k + 2.
942 S. TAŞPINAR ET AL.

found in Anselin (1988), Pinkse and Slade (1998), Fingleton (2009), Burridge and Fingleton (2010), and
Jin and fei Lee (2013, 2015). According to the bootstrap approach, the empirical distribution of a test
statistic is constructed by resampling data. For our model, sampling from Yn , Wn Yn , and Xn cannot
be done without preserving the inherent spatial dependence structure of the model. Instead, Anselin
(1988) considers bootstrapping of residuals when disturbances  of the model are i.i.d. Given a consistent
estimator θ̂n , the vector of residuals is ε̂n = Rn (ρ̂n ) Sn (λ̂n )Yn −Xn β̂n . Under Assumption 1, we demean
residuals to obtain ε̃n = Jn ε̂n , where Jn = In − n1 ln l′ n . The bootstrap sample of residuals, denoted by
εn∗ , is generated with random draws from the empirical distribution of the demeaned residuals. Then,
 
the vector of random draws εn∗ is used to generate Yn∗ = S−1 −1 ∗
n (λ̂n ) Xn β̂n + R (ρ̂n )εn . As a result, a
∗,X ′ ′
bootstrap sample of (Yin n,i• ) : i = 1, . . . , n , where Xn,i• is the ith row of Xn , is obtained. For our
GMM set up, the bootstrapped Wald test of H0 : λ0 = 0 is implemented with the following  procedure:
(i) Compute the best GMME θ̂n stated in (3.3), and the resulting ε̃n = Jn Rn (ρ̂n ) Sn (λ̂n )Yn − Xn β̂n .
Calculate tλ = λ̂n /Se(λ̂n ), where Se(λ̂n ) is the estimated asymptotic standard errors of λ̂ n .
(ii) Draw εn⋆ from  ε̃n through sampling with replacement. Generate Yn∗ = S−1 n (λ̂n ) Xn β̂n +
−1
R (ρ̂n )εn . ∗

(iii) Formulate an initial GMME with the bootstrap version of moment functions: gn∗ (θ ) =
′ ′ ′


εn∗ (θ )P1n εn∗ (θ ), . . . , εn∗ (θ )P4n εn∗ (θ ), εn∗ (θ )Qn ′ , where P1n = Wn ′ Wn − D Wn ′ Wn , P2n =


(t)
(t)

Mn ′ Mn −D Mn ′ Mn , P3n = Wn ′ Wn , P4n = Mn ′ Mn , Qn = Wn Mn Xn , Wn Xn , Mn Xn , Xn ,
′
and εn∗ (θ ) = Rn (ρ)Sn (λ)Yn∗ − Rn (ρ)Xn β. Define θ̂b1 = argminθ gn∗ (θ )gn∗ (θ ).
(iv) Using θ̂b1 and corresponding residuals, construct an estimate of ⋆n given in (3.1) and denote it
with  ˆ ∗n . Using θ̂b1 , generate the best set of moment functions gn∗ (θ ). Define the bootstrap optimal
′
GMME by θ̂b2 = argminθ gn∗ (θ )  ˆ ∗−1 ∗
n gn (θ ).
(v) Using θ̂b2 and corresponding residuals, find an estimate of 8⋆n and denote it with 8̂∗n . Obtain an
′ ∗−1 ∗ −1
estimate of the variance of θ̂b2 from (8̂∗n  ˆ n 8̂n ) . Obtain the bootstrap standard error of λ̂b2


and denote it with Se(λ̂b2 ). Then, calculate tλb = λ̂2b − λ̂n /Se(λ̂b2 ).
P

(vi) Repeat steps (ii)–(v) B times and calculate the bootstrap p-value of pλ = B1 Bb=1 1 |tλb | > |tλ | .
Then, the test rejects the null at nominal size α if pλ < α.
As for a test of H0 : ρ0 = 0, simply compute tρ in step (i) and tρb in step (v), and evaluate these values
similar to step (vi).
Given the result in (3.4), tλ has an asymptotic standard normal distribution, and hence the inference
is based on this asymptotic distribution. On the other hand, the distribution function of tλb obtained from
the bootstrap method is an approximation to the exact finite sample distribution function of tλ and can
differ from the asymptotic distribution.13

6. A Monte Carlo study

6.1. Design
For the model
 specified in (2.1), we assume that there are are two regressors and no intercept term such
that Xn = X1,n , X2,n and β0 = (β10 , β20 )′ = (1, 1)′ , where X1,n and X2,n are n × 1 vectors of variables.
The i-th observation of disturbance term εin is i.i.d. Normal(0, 1). The spatial autoregressive parameters
λ0 and ρ0 take values from the set K = {−0.6, −0.3, 0, 0.3, 0.6} in order to allow for both weak and
strong spatial interactions.
As for X1,n and X2,n , we use the county-level data set of Pace and Barry (1997) on the 1980 presidential
election. We set X1,n equal to the standardized value of income per-capita and X2,n equal to the

13
The Wald test is not the only test statistic that could be used for hypothesis testing. We refer the reader to Web Appendix
C for the formulation of five other test statistics (the LM test, the C(α) test, the minimum chi-square test (MC), and two
criterion-based tests (DNW , and DRU )), and we report on their performance in Web Appendix E.
 ECONOMETRIC REVIEWS 943

standardized value of the homeownership rate. The data set describes 3,107 U.S. counties, of which we
use the first n observations in the Monte Carlo study.
The weight matrix Wn is based on the interaction scenario described in Arraiz et al. (2010).
The observations are distributed across four quadrants of a space in such a way that the number of
observations in each quadrant can be arranged to allow for sparse or dense quadrants. The location
of each unit across the space is determined by the xy-coordinates over a square grid. Let m and m̄ be
two integers. Then, the units in the northeast quadrant of the space have discrete coordinates satisfying
m + 1 ≤ x ≤ m̄ and m + 1 ≤ y ≤ m̄, with an increment value of 0.5. For the other quadrants, the
location coordinates are integers satisfying 1 ≤ x ≤ m, 1 ≤ y ≤ m̄, and 1 ≤ x ≤ m̄, 1 ≤ y ≤ m.
The distance dij between any two units i and j, located respectively at (x1 , y1 ) and (x2 , y2 ), is measured
 1/2
by the Euclidean distance given by dij = (x1 − x2 )2 + (y1 − y2 )2 . Then, the (i, j)th element of the
Wn , wn,ij , equals 1 if 0 ≤ dij ≤ 1 and zero otherwise.
Varying the values for m and m̄ leads to a different sample size and a different share of units in
the dense northeast quadrant. For our small samples, we consider combinations (m, m̄) = (5, 15)
and (m, m̄) = (14, 20). The first combination produces a sample size of 486 and locates 75% of
units in the northeast quadrant, whereas the second combination generates a sample size of 485 and
locates 25% of the observation in the northeast quadrant. To generate large sample sizes, we specify
(m, m̄) = (7, 21) and (m, m̄) = (20, 28). Here, the first combination generates a sample size of 974
with 75% of observation located in the northeast quadrant, and the second one produces a sample
size of 945 with 25% units in the northeast quadrant. Let the percentage of units located in the
northeast quadrant be denoted by d. Then our Monte Carlo experiments are based on the combinations
(m, m̄, n, d) = {(5, 15, 486, 75%), (14, 20, 485, 25%), (7, 21, 974, 75%), (20, 28, 945, 25%)}, and we refer
to these combinations by their (n, d) values.
For the spatial weight matrix Mn , we consider a distance-based binary matrix. Let Ji be the set of the
nearest 10 observations to the ith observation based on dij . Then, the (i, j)th element of Mn , mn,ij , equals
1 if j ∈ Ji and zero otherwise.
Both Wn and Mn are row normalized. For each specification, resampling is repeated 1,000 times. We
examine the following estimators: (i) the IGMME, (ii) the OGMME, and (iii) the MLE. For each esti-
mator, we provide bias, estimated asymptotic standard errors (OGMME-Asy.), and empirical standard
deviations (Emp.). Besides these results, we also provide the corrected standard errors (OGMME-Crc.)
in the case of OGMME.14

6.2. Estimation results

The results are easily interpretable when presented in graphical form, and in doing so, we highlight the
standard error estimators. To this end, we compute the percentage deviation of the estimated asymptotic
and corrected standard errors from the corresponding empirical standard deviations; i.e., we compute
100 × (Asy. − Emp.)/Emp. and 100 × (Crc. − Emp.)/Emp. We then plot these percentage deviations
against 25 combinations of (λ0 , ρ0 ) separately for each (n, d) value. Following Bond and Windmeijer
(2005, p. 18), we use the term “bias” for these deviations. If the asymptotic distribution of an estimator
performs well enough, the bias should be near zero. Conversely, a large bias suggests that, for the given
estimator, the asymptotic distribution approximates the finite sample distribution poorly.
Figure 1 presents the results for ρ̂n . There is a downward bias in the case of OGMME-Asy., OGMME-
Crc., and MLE. The direction of the bias for the IGMME changes from positive to negative as the
northeast quadrant density changes from low to high. The estimated asymptotic standard errors of
IGMME can exceed the corresponding empirical standard errors by more than 10%. In the case of the
MLE, spikes for ρ0 = −0.6 result in a downward bias of up to 50%. Except for these spikes, the estimated
asymptotic standard errors of the MLE are relatively closer than those of OGMME to the corresponding

14
Asy. and Crc. are the median values of their respective distributions. For the MLE, we use the analytical result for the
covariance matrix to calculate standard errors.
944 S. TAŞPINAR ET AL.

Figure 1. Percentage deviations of the estimated standard errors from the empirical standard deviations: ρ̂n .

empirical standard deviations. In the case of OGMME, the correction method provides standard errors
that are closer to the empirical standard deviations. In Fig. 1(b), where the density of the northeast
quadrant is high, all estimators have a downward bias. The size of the bias is the smallest in the case of
the correction method for all combinations of (λ0 , ρ0 ). In Figs. 1(c) and 1(d), where the results for the
large samples are presented, the same pattern is also observed. That is, the standard errors based on the
correction method are relatively closer to the corresponding empirical standard deviations.
The results for λ̂n are presented in Fig. 2. In Figs. 2(a) and 2(b), the magnitude of the percentage
deviations are similar for the OGMME-Asy. and the OGMME-Crc. In both figures, the bias is generally
less than 5%, and it is relatively larger in Fig. 2(a). As can be seen from these figures, the correction
method generally provides an improvement for the statistical inference as they are closer to the empirical
standard deviations. Again for the combinations of autoregressive parameters where ρ0 = −0.6, the
MLE results display spikes, but this time with smaller amplitudes than for ρ̂n . The asymptotic method
based on the IGMME provides standard error estimates that are biased by 25% or more. As the northeast
quadrant density increases, biases become smaller. Biases are generally less than 5% in the case of
OGMME-Asy., OGMME-Crc., and MLE for all (λ0 , ρ0 ) combinations.
 ECONOMETRIC REVIEWS 945

Figure 2. Percentage deviations of the estimated standard errors from the empirical standard deviations: λ̂n .

We summarize the main findings of our simulation experiment in the following list with a focus on
results reported for the OGMME-Crc.15 :
(i) Across the various (λ0 , ρ0 ) combinations, the three estimators of θ0 have a negligible bias.
(ii) The size of the northeast quadrant in our experiment affects the precision of estimators. When d
is larger, the empirical standard deviations of the IGMME for ρ̂n and λ̂n are relatively larger. The
MLE for λ̂n , β̂1n , and β̂2n has relatively smaller empirical standard deviations, and the standard
errors of the OGMME-Crc. for ρ̂n is relatively smaller.
(iii) For ρ̂n , when n = 485 or n = 486, our results indicate that the downward bias in the corrected
standard errors is less than 5%, while the downward bias in the estimated asymptotic standard
errors is generally larger than 5%. When n = 945 or n = 974, the corrected standard errors still
perform better than the estimated asymptotic standard errors and again they are less than 5% in
almost all cases. In order to give an overall picture, we compute for each estimator the proportions
of percentage deviations that fall within a (−5%, 5%) interval across all combinations of (λ0 , ρ0 )

15
The results for all cases are provided in tables presented in Web Appendix D.
946 S. TAŞPINAR ET AL.

and (n, d) values. These proportions are (i) 0.33 for the IGMME, (ii) 0.40 for the OGMME-Asy.,
(iii) 0.89 for the OGMME-Crc., and (iv) 0.67 for the MLE.
(iv) The corrected standard errors for λ̂n has relatively less downward bias than the estimated
asymptotic standard errors. If we consider the proportions of percentage deviations that fall
within a (−5%, 5%) interval for each estimator, we have (i) 0.27 for the IGMME, (ii) 0.80 for
the OGMME-Asy., (iii) 0.88 for the OGMME-Crc., and (iv) 0.84 for the MLE.
(v) For β̂1n , both the OGMME-Asy. and the OGMME-Crc. have a downward bias of more than 5% in
some cases. The proportion of percentage deviations that fall within a (−5%, 5%) interval equals
(i) 0.57 for the IGMME, (ii) 0.36 for the OGMME-Asy., (iii) 0.24 for the OGMME-Crc., and (iv)
0.93 for the MLE. It is obvious that the MLE is performing best in this case.
(vi) On β̂2n , both the OGMME-Asy. and OGMME-Crc. provide similar results with a bias amount less
than 5% in many cases. The proportion of percentage deviations that fall in a (−5%, 5%) interval
is (i) 0.48 for the IGMME, (ii) 0.83 for the OGMME-Asy., (iii) 0.95 for the OGMME-Crc., and
(iv) 0.92 for the MLE.

6.3. Inference results

For each test, we consider the empirical rejection frequencies under the null hypothesis of the form H0 :
λ0 = 0 or H0 : ρ0 = 0. We use p-value plots suggested by Davidson and MacKinnon (1998) to illustrate
the actual size of these tests under the null hypothesis for various nominal size values. Figures 3 and 4
are the p-value plots for the null hypothesis of no spatial dependence, respectively, in the disturbance
term and in the dependent variable. We present the results for the Wald tests and its bootstrap version
for two representative cases and leave the results for other cases and other test statistics to Web Appendix
E. In the figures, W1 , W2 , W2c , Wm , and W2b denote respectively the Wald tests based on the IGMME,
OGMME-Asy, OGMME-Crc, MLE, and the bootstrap version.
Conventional significance levels in hypothesis tests are often chosen to be 10% or smaller. Therefore,
we use nominal size values that are less than or equal to 12% in these figures. If the asymptotic
approximation used for the distributions of these test statistics perform well, we would expect that the
line graphs in these figures to follow closely the 45◦ line. Hence, a discrepancy from the 45◦ line indicates
that the sampling distribution of the test statistic differs from the asymptotic distribution that is used to
establish the critical values.
First, we evaluate properties of tests under the null of H0 : ρ0 = 0. For the null hypothesis under
consideration, we only present the p-value plots for the case of λ0 = 0.3 in Fig. 3. The p-value plots for
all other cases are provided in the web appendix. The salient properties of these tests are summarized as
follows:
(i) The Wald test based on the IGMME (W1 ) in Fig. 3 is not properly sized. It is undersized when
the northeast density is low and is over-sized when the northeast density is high. The increase in
the sample size does not improve the size property of W1 but rather makes it worse for the low
density case.
(ii) The standard Wald test based on the OGMME-Asy (W2 ) is over-sized in all cases which confirms
previous findings in the literature on non-spatial econometric models. Therefore, inference based
on this test could be very misleading. For example, at the nominal size of 6%, this test falsely reject
the null in around 10% of the replications. This test performs worse than the Wm and W2c tests in
all cases. This result is principally due to the downward bias in the estimated asymptotic standard
errors of the OGMME reported in Section 6.2.
(iii) For large sample sizes of n = 945 and n = 974, W2 seems to improve but is over-sized relative
to Wm and W2c tests. The size of the northeast quadrant appears to affect the performance of this
test in small sample sizes.
(iv) The corrected version based on the OGMME-Crc (W2c ) outperforms W2 in all cases, which
confirms our results reported in Section 6.2. Moreover, it performs as well as the MLE-based
 ECONOMETRIC REVIEWS 947

Figure 3. Wald tests: H0 : ρ0 = 0 (for λ0 = 0.3).

test (Wm ). Thus, the correction method improves the size properties of the Wald test, which is
consistent with the simulation results presented in Section 6.2.
(v) The MLE based Wald test (Wm ) has a proper size in all cases: the asymptotic distribution derived
under the null hypothesis provides a good approximation to the finite sample distribution.
(vi) As the bootstrapped version of the Wald test (W2b ) is computationally intensive, we have
computed results only for the OGMME and only for the cases where n = 485 and n = 486.16 The
W2b test performs better than the W2 test and is slightly over-sized for larger nominal size values.
In comparison with the Wm and W2c tests, the size is slightly worse.
(vii) The size properties of remaining tests, LM, MC, DNW , DRU , and C(α), are evaluated in detail in
Web Appendix.17 Overall, the DNW and DRU tests perform better than the other tests. The C(α)
test is properly sized in many cases and performs better than LM and MC.
Next, we evaluate the size properties of these tests under the null of H0 : λ0 = 0. For the null
hypothesis under consideration, we only present the p-value plots for the case of ρ0 = 0.3 in Fig. 4.

16
For each resample, the number of bootstrap samples is limited to 99. This choice for the number of bootstrap samples
satisfies the rule of thumb provided by Davidson and MacKinnon (1998): (99 + 1) × 0.01 is a positive integer. A larger
number of bootstrap samples is computationally challenging for the spatial model under consideration.
17
Here, MC denotes the minimum chi-square test, and DNW and DRU denote criterion based tests. For details, see Web
Appendix C.
948 S. TAŞPINAR ET AL.

Figure 4. Wald tests: H0 : λ0 = 0 (for ρ0 = 0.3).

The p-value plots for all other cases are provided in the Web Appendix. The salient properties of these
tests are as follows:
(i) The W1 test is oversized in all cases regardless of the density of the northeast quadrant and the
sample size. For example, in Fig. 4(d), when the nominal size is 6% this test rejects the true null
more than 10% of the replications. As shown in Figs. 4(a) to 4(d), the rejection frequencies do not
converge to the true ones as the sample size increases.
(ii) The standard Wald test W2 is slightly oversized at low nominal sizes when d = 0.25 and more
generally so when d = 0.75.
(iii) The corrected version W2c achieves the correct sizes implied by the asymptotic theory in small
samples. In large samples, it is slightly undersized when d = 0.25 and slightly oversized when
d = 0.75.
(iv) The results for the bootstrapped version are again provided only for the Wald test with the
OGMME when n = 485 and n = 486. As seen, this test is over-sized and its deviation from
the 45◦ line becomes larger for the large nominal size values. For example, for the nominal sizes
of 6% and 10% the test statistic produces actual sizes of 9% and 11%, respectively in Fig. 4(a).
Note, though, that the number of bootstrap samples is limited to only 99: this evidence may yet
be subject to random fluctuations in the resampling process.
 ECONOMETRIC REVIEWS 949

(v) The MLE based Wald test is undersized when (n, d) = (485, 25%) and is slightly over-sized when
(n, d) = (945, 25%) for larger nominal values. As the northeast quadrant density increases, the
actual sizes converge to the nominal sizes.
(vi) The size properties of remaining tests are examined in the Web Appendix. Both LM and MC
tests are over-sized in all cases and are not recommended for applied research to test this null
hypothesis. The two criterion based tests, DNW and DRU , are properly sized only when there exist
negative spatial dependence in the disturbance term. The p-value plots for the C(α) test are close
to the 45◦ line in many cases, indicating that this test is properly sized.

6.4. Results on the effects estimators

We use the same Monte Carlo setup described in Section 6.1 to study the finite sample properties
of the estimators of the impact measures defined in Section 4. For the sake of brevity, we consider
only (n, d) = {(486, 75%), (485, 46%), (485, 25%), (974, 75%), (945, 46%), (945, 25%)}. Besides the
IGGME, OGMME, and MLE, we also consider the performance of the multi-step estimator (GS2SLSE)
of Drukker et al. (2013) in this analysis.18 For each impact measure, we provide the bias, the estimated

Table 1. Average absolute bias of impact measures.

Total Direct Indirect
Simulation Delta Simulation Delta Simulation Delta
X1,n X2,n X1,n X2,n X1,n X2,n X1,n X2,n X1,n X2,n X1,n X2,n
n = 486, d = 75%
IGMME 0.0194 0.0171 0.0091 0.0080 0.0029 0.0032 0.0024 0.0028 0.0192 0.0163 0.0088 0.0070
OGMME_Asy. 0.0213 0.0183 0.0090 0.0091 0.0034 0.0035 0.0031 0.0036 0.0221 0.0173 0.0098 0.0075
MLE 0.0160 0.0131 0.0099 0.0092 0.0022 0.0030 0.0017 0.0031 0.0142 0.0117 0.0093 0.0079
GS2SLSE 0.0378 0.0331 0.0215 0.0195 0.0035 0.0034 0.0025 0.0031 0.0355 0.0306 0.0206 0.0181
n = 485, d = 46%
IGMME 0.0112 0.0075 0.0064 0.0095 0.0028 0.0050 0.0028 0.0040 0.0104 0.0082 0.0064 0.0079
OGMME_Asy. 0.0145 0.0108 0.0060 0.0061 0.0028 0.0037 0.0030 0.0036 0.0137 0.0110 0.0048 0.0061
MLE 0.0087 0.0077 0.0062 0.0059 0.0024 0.0033 0.0027 0.0030 0.0066 0.0060 0.0070 0.0056
GS2SLSE 0.0276 0.0216 0.0143 0.0114 0.0034 0.0046 0.0024 0.0034 0.0258 0.0215 0.0143 0.0118
n = 485, d = 25%
IGMME 0.0140 0.0084 0.0075 0.0071 0.0039 0.0054 0.0033 0.0044 0.0146 0.0094 0.0077 0.0071
OGMME_Asy. 0.0182 0.0127 0.0101 0.0065 0.0034 0.0032 0.0030 0.0030 0.0167 0.0128 0.0078 0.0062
MLE 0.0130 0.0096 0.0053 0.0059 0.0030 0.0030 0.0027 0.0028 0.0131 0.0093 0.0062 0.0053
GS2SLSE 0.0335 0.0253 0.0179 0.0141 0.0056 0.0054 0.0035 0.0039 0.0317 0.0266 0.0182 0.0160
n = 974, d = 75%
IGMME 0.0112 0.0063 0.0056 0.0048 0.0030 0.0035 0.0029 0.0031 0.0108 0.0075 0.0057 0.0050
OGMME_Asy. 0.0099 0.0052 0.0063 0.0050 0.0023 0.0028 0.0024 0.0028 0.0094 0.0063 0.0049 0.0045
MLE 0.0059 0.0042 0.0062 0.0061 0.0025 0.0023 0.0024 0.0025 0.0049 0.0038 0.0061 0.0055
GS2SLSE 0.0164 0.0125 0.0082 0.0071 0.0029 0.0035 0.0024 0.0030 0.0172 0.0134 0.0091 0.0080
n = 945, d = 47%
IGMME 0.0087 0.0048 0.0077 0.0055 0.0031 0.0030 0.0026 0.0028 0.0098 0.0055 0.0075 0.0062
OGMME_Asy. 0.0083 0.0052 0.0056 0.0038 0.0024 0.0030 0.0020 0.0027 0.0091 0.0061 0.0053 0.0038
MLE 0.0055 0.0046 0.0044 0.0046 0.0022 0.0026 0.0017 0.0024 0.0053 0.0055 0.0049 0.0047
GS2SLSE 0.0134 0.0106 0.0060 0.0054 0.0036 0.0031 0.0026 0.0025 0.0131 0.0111 0.0066 0.0061
n = 945, d = 25%
IGMME 0.0115 0.0058 0.0048 0.0049 0.0026 0.0030 0.0020 0.0024 0.0102 0.0053 0.0054 0.0042
OGMME_Asy. 0.0079 0.0049 0.0038 0.0028 0.0020 0.0024 0.0021 0.0024 0.0079 0.0060 0.0038 0.0033
MLE 0.0064 0.0040 0.0041 0.0032 0.0019 0.0023 0.0018 0.0024 0.0052 0.0042 0.0040 0.0037
GS2SLSE 0.0112 0.0089 0.0052 0.0045 0.0029 0.0029 0.0020 0.0023 0.0098 0.0081 0.0045 0.0039

18
Drukker et al. (2013) extend the GMM/IV estimation method of Kelejian and Prucha (1998) to a SARAR(1,1) specification
that includes endogenous regressors. Although they study the finite sample properties of estimators and the Wald tests,
the performance of impact measures is not explored. For our discussion of their estimators, see Web Appendix F.
950 S. TAŞPINAR ET AL.

asymptotic standard errors, the empirical standard deviations, and the rejection rates of Wald test for a
nominal size value of 0.05.
The simulation results are given in 36 tables presented in Web Appendix G. Here, we summarize these
results by providing scalar measures. We start by comparing the bias properties of effect estimates across
estimators. To evaluate the performance of estimators in terms of bias, we calculate the average absolute
bias for each impact measure over 25 combinations of (λ0 , ρ0 ). These results are given in Table 1. In
general, the average absolute bias reported by each estimator is less than 0.05 for each impact measure.
The IGMME and GS2SLSE impose slightly larger bias on impact measures. Table 1 also indicates that
the Delta method performs slightly better than the simulation method. The magnitude of bias is slightly
larger when the density of the northeast quadrant is high.
Next, we compare the estimated standard errors with the empirical standard deviations for each
estimator. For the total effects, we provide the percentage deviations of the estimated standard errors
from the empirical standard deviations in Fig. 5 for the case of (n, d) = (486, 75%). Generally, the
percentage bias falls within the (−5%, 5%) interval for each estimator, except in the case of IGMME. To
give an overall picture, we calculate the proportions of percentage deviations that fall within (−5%, 5%)

Figure 5. Total effects: Percentage deviations of the estimated standard errors from the empirical standard deviations for n = 486,
d = 75%.
 ECONOMETRIC REVIEWS 951

across all combinations of (λ0 , ρ0 ) and (n, d). For the simulation method, these proportions for the total
effects of X1,n are 0.48 for the IGMME, 0.746 for the OGMME-Asym., 0.68 for the OGMME-Crc., 0.94
for the MLE, and 0.933 for the GS2SLSE, and for X2,n these rates are 0.526 for the IGMME, 0.80 for
the OGMME-Asym., 0.88 for the OGMME-Crc., 0.953 for the MLE, and 0.933 for the GS2SLSE. For
the Delta method, these proportions for the total effects of X1n are 0.453 for the IGMME, 0.70 for the
OGMME-Asym., 0.74 for the OGMME-Crc., 0.913 for the MLE, and 0.926 for the GS2SLSE, and for X2n
these rates are 0.506 for the IGMME, 0.873 for the OGMME-Asym., 0.933 for the OGMME-Crc., 0.92
for the MLE, and 0.94 for the GS2SLSE. While the IGMME performs relatively poorly, the MLE and
the GS2SLSE perform relatively better. In general, the OGMME-Crc performs better than OGMME-
Asy, which is consistent with our results in Sections 6.2 and 6.3. The simulation method seems to be
performing slightly better than the Delta method in many cases. The standard errors of impact estimators
seem to be affected by the density of northeast quadrant, and the proportions of percentage deviations
that falls in the (−5%, 5%) interval are relatively higher for the case of (n, d) = (485, 46%).

Table 2. Average size distortions of the Wald tests for impact measures.
Total Direct Indirect
Simulation Delta Simulation Delta Simulation Delta
X1,n X2,n X1,n X2,n X1,n X2,n X1,n X2,n X1,n X2,n X1,n X2,n
n = 486, d = 75%
IGMME 0.0187 0.0216 0.0224 0.0266 0.0107 0.0124 0.0107 0.0121 0.0272 0.0249 0.0314 0.0294
OGMME_Asy. 0.0132 0.0092 0.0165 0.0112 0.0210 0.0103 0.0217 0.0107 0.0110 0.0106 0.0137 0.0142
OGMME_Crc. 0.0065 0.0061 0.0068 0.0072 0.0252 0.0069 0.0280 0.0073 0.0065 0.0102 0.0090 0.0120
MLE 0.0082 0.0085 0.0100 0.0118 0.0070 0.0063 0.0070 0.0063 0.0080 0.0087 0.0106 0.0106
GS2SLSE 0.0054 0.0059 0.0064 0.0070 0.0068 0.0058 0.0064 0.0061 0.0064 0.0057 0.0059 0.0066
n = 485, d = 46%
IGMME 0.0136 0.0138 0.0140 0.0154 0.0102 0.0160 0.0100 0.0153 0.0244 0.0220 0.0264 0.0252
OGMME_Asy. 0.0127 0.0099 0.0146 0.0106 0.0192 0.0115 0.0193 0.0112 0.0087 0.0086 0.0095 0.0105
OGMME_Crc. 0.0066 0.0073 0.0058 0.0075 0.0212 0.0073 0.0232 0.0074 0.0050 0.0081 0.0056 0.0075
MLE 0.0074 0.0073 0.0096 0.0076 0.0080 0.0098 0.0084 0.0089 0.0065 0.0079 0.0079 0.0095
GS2SLSE 0.0066 0.0074 0.0070 0.0085 0.0076 0.0085 0.0078 0.0077 0.0051 0.0055 0.0055 0.0068
n = 485, d = 25%
IGMME 0.0222 0.0159 0.0233 0.0192 0.0116 0.0153 0.0123 0.0154 0.0392 0.0342 0.0415 0.0387
OGMME_Asy. 0.0118 0.0078 0.0143 0.0093 0.0186 0.0083 0.0180 0.0088 0.0083 0.0070 0.0106 0.0091
OGMME_Crc. 0.0079 0.0051 0.0076 0.0066 0.0274 0.0062 0.0282 0.0067 0.0048 0.0075 0.0059 0.0082
MLE 0.0109 0.0092 0.0122 0.0114 0.0078 0.0080 0.0079 0.0081 0.0105 0.0102 0.0126 0.0122
GS2SLSE 0.0060 0.0063 0.0076 0.0076 0.0058 0.0050 0.0051 0.0048 0.0079 0.0065 0.0078 0.0068
n = 974, d = 75%
IGMME 0.0212 0.0176 0.0218 0.0176 0.0159 0.0218 0.0157 0.0220 0.0278 0.0223 0.0273 0.0223
OGMME_Asy. 0.0103 0.0069 0.0108 0.0076 0.0112 0.0084 0.0109 0.0086 0.0092 0.0068 0.0099 0.0078
OGMME_Crc. 0.0056 0.0062 0.0058 0.0062 0.0128 0.0071 0.0135 0.0070 0.0077 0.0065 0.0075 0.0071
MLE 0.0065 0.0052 0.0069 0.0063 0.0032 0.0051 0.0036 0.0050 0.0063 0.0070 0.0074 0.0078
GS2SLSE 0.0045 0.0051 0.0048 0.0058 0.0036 0.0064 0.0037 0.0059 0.0058 0.0049 0.0044 0.0048
n = 945, d = 47%
IGMME 0.0265 0.0160 0.0268 0.0163 0.0144 0.0191 0.0148 0.0194 0.0406 0.0344 0.0419 0.0361
OGMME_Asy. 0.0137 0.0091 0.0137 0.0090 0.0132 0.0066 0.0129 0.0072 0.0080 0.0084 0.0096 0.0093
OGMME_Crc. 0.0074 0.0059 0.0075 0.0066 0.0150 0.0063 0.0147 0.0064 0.0076 0.0080 0.0070 0.0083
MLE 0.0068 0.0074 0.0083 0.0078 0.0064 0.0044 0.0058 0.0043 0.0056 0.0049 0.0064 0.0054
GS2SLSE 0.0064 0.0069 0.0067 0.0077 0.0056 0.0039 0.0058 0.0041 0.0059 0.0061 0.0061 0.0062
n = 945, d = 25%
IGMME 0.0194 0.0156 0.0189 0.0152 0.0153 0.0220 0.0156 0.0222 0.0319 0.0264 0.0322 0.0263
OGMME_Asy. 0.0094 0.0082 0.0100 0.0080 0.0122 0.0070 0.0115 0.0066 0.0069 0.0056 0.0071 0.0059
OGMME_Crc. 0.0050 0.0070 0.0049 0.0066 0.0137 0.0054 0.0127 0.0051 0.0054 0.0052 0.0055 0.0056
MLE 0.0054 0.0060 0.0055 0.0065 0.0059 0.0065 0.0057 0.0054 0.0051 0.0062 0.0056 0.0065
GS2SLSE 0.0042 0.0045 0.0042 0.0044 0.0052 0.0047 0.0054 0.0046 0.0038 0.0038 0.0039 0.0034
952 S. TAŞPINAR ET AL.

In Table 2, we report the average size distortions of the Wald tests over 25 combinations of (λ0 , ρ0 )
for each impact measure for a nominal size of 0.05.19
The average size distortions are smaller than 0.03 in all cases, and the Wald tests based on IGMME
have relatively larger average size distortions. Overall, the average size distortions of the Wald tests based
on OGMME-Crc are smaller than those based on OGMME-Asy, which suggests that the correction
method can lead to more accurate inference.

7. Conclusion
In this study, we review the GMM estimation approach for spatial autoregressive models and evaluate
the finite sample properties of the GMM estimators through simulation studies. To the best of our
knowledge, this study is the first extensive study to evaluate various GMM inference techniques for
spatial autoregressive models. We extend the finite sample correction method of Windmeijer (2005) to
the SARAR(1,1) model and formulate a variance formula to improve inference properties of the GMM
estimator for detecting the presence of spatial dependence. We show that the corrected standard errors
can be a good alternative to their estimated asymptotic counterparts.
Simulation studies confirm that the estimated asymptotic standard errors based on the OGMME can
be substantially downward biased and can make inference problematic for the spatial autoregressive
parameters. This result is consistent with the simulation results reported in the econometrics literature
on the nonspatial models. More specifically, we evaluate the percentage deviations of the corrected
standard errors from the empirical standard deviations under various scenarios. The results show
that the bias in the corrected standard errors of OGMME and the effects estimators (or the impacts
estimators) is mostly less than the bias in the estimated asymptotic standard errors of OGMME and
the effects estimators. Hence, the correction method can provide more accurate inference in finite
samples.
For the SARAR(1,1) model, we also compare the size properties of various test statistics through
simulations. These tests included standard asymptotic Wald tests based on MLE and GMMEs, a
bootstrapped version of Wald test, two versions of C(α) test, the standard LM test, the minimum chi-
square test, and two versions of GMM criterion test. Our results show that the standard asymptotic Wald
test based on the asymptotic standard errors of the OGMME is generally oversized, whereas the Wald test
based on the corrected standard errors of the OGMME is properly sized. Our findings also confirm that
the GMM criterion tests and the C(α) test can be useful for detecting the presence spatial dependence.
All of these findings offer guidance to applied researchers who estimate and test spatial models with the
GMM estimators.

Acknowledgments
We thank the editor, an associate editor, two anonymous referees and the conference participants at SEA2015 (Miami) and
Camp Econometrics X (New York).

Funding
This research was supported, in part, by a grant of computer time from the City University of New York High Performance
Computing Center under NSF Grants CNS-0855217 and CNS-0958379.

19
We calculate the average size distortions over 25 combinations of (λ0 , ρ0 ) from 36 tables given in Web Appendix G by
1 P25 |actual size − 0.05|, where i corresponds to the ith combinations of (λ , ρ ).
25 i=1 i 0 0
 ECONOMETRIC REVIEWS 953

References
Altonji, J. G., Segal, L. M. (1996). Small-sample bias in gmm estimation of covariance structures. Journal of Business &
Economic Statistics 14(3):353–366.
Anselin, L. (1988). Spatial econometrics: Methods and Models. New York: Springer.
Arellano, M., Bond, S. (1991). Some tests of specification for panel data: Monte carlo evidence and an application to
employment equations. The Review of Economic Studies 58(2):277–297.
Arraiz, I. et al. (2010). A spatial Cliff-Ord-type model with heteroskedastic innovations: Small and large sample results.
Journal of Regional Science 50(2):592–614.
Blundell, R., Bond, S. (1998). Initial conditions and moment restrictions in dynamic panel data models. Journal of
Econometrics 87(1):115–143.
Bond, S., Windmeijer, F. (2005). Reliable inference for gmm estimators? finite sample properties of alternative test
procedures in linear panel data models. Econometric Reviews 24(1):1–37.
Brown, B. W., Newey, W. K. (2002). Generalized method of moments, efficient bootstrapping, and improved inference.
Journal of Business & Economic Statistics 20(4):507–517.
Burridge, P., Fingleton, B. (2010). Bootstrap inference in spatial econometrics: The J-test. Spatial Economic Analysis
5(1):93–119.
Clark, T. E. (1996). Small-sample properties of estimators of nonlinear models of covariance structure. Journal of Business
& Economic Statistics 14(3):367–373.
Das, D., Kelejian, H. H., Prucha, I. R. (2003). Small sample properties of estimators of spatial autoregressive models with
autoregressive disturbances. Papers in Regional Science 82:1–26.
Davidson, R., MacKinnon, J. G. (1998). Graphical methods for investigating the size and power of hypothesis tests. The
Manchester School 66(1):1–26.
Debarsy, N., Jin, F., Lee, L. f. (2015). Large sample properties of the matrix exponential spatial specification with an
application to FDI. Journal of Econometrics 188(1):1–21.
Drukker, D. M., Egger, P., Prucha, I. R. (2013). On two-step estimation of a spatial autoregressive model with autoregressive
disturbances and endogenous regressors. Econometric Reviews 32(5–6):686–733.
Fingleton, B. (2009). A generalized method of moments estimator for a spatial model with moving average errors, with
application to real estate prices. In: Arbia, G., Baltagi, B. H., eds. Spatial Econometrics. Studies in Empirical Economics,
pp. 35–57.
Hall, P., Horowitz, J. L. (1996). Bootstrap critical values for tests based on generalized-method-of-moments estimators.
Econometrica 64(4):891–916.
Hansen, L. P., Heaton, J., Yaron, A. (1996). Finite-sample properties of some alternative gmm estimators. Journal of Business
& Economic Statistics 14(3):262–280.
Jin, F., Lee, L. f. (2013). Cox-type tests for competing spatial autoregressive models with spatial autoregressive disturbances.
Regional Science and Urban Economics 43(4):590–616.
Jin, F., Lee, L. f. (2015). On the bootstrap for moran‘s test for spatial dependence. Journal of Econometrics 184(2):295–314.
Kelejian, H. H., Prucha, I. R. (1998). A generalized spatial two-stage least squares procedure for estimating a spatial
autoregressive model with autoregressive disturbances. Journal of Real Estate Finance and Economics 17(1):1899–1926.
Kelejian, H. H., Prucha, I. R. (2010). Specification and estimation of spatial autoregressive models with autoregressive and
heteroskedastic disturbances. Journal of Econometrics 157:53–67.
Lee, L.-f. (2003). Best spatial two-stage least squares estimators for a spatial autoregressive model with autoregressive
disturbances. Econometric Reviews 22(4):307–335.
Lee, L.-f. (2004). Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models.
Econometrica 72(6):1899–1925.
Lee, L.-f. (2007a). GMM and 2SLS estimation of mixed regressive, spatial autoregressive models. Journal of Econometrics
137(2):489–514.
Lee, L.-f. (2007b). The method of elimination and substitution in the GMM estimation of mixed regressive, spatial
autoregressive models. Journal of Econometrics 140:155–189.
Lee, L.-f., Liu, X. (2010). Efficient GMM estimation of high order spatial autoregressive models with autoregressive
disturbances. Econometric Theory, 26(1):187–230.
LeSage, J., Pace, R. K. (2009). Introduction to Spatial Econometrics. Statistics: A Series of Textbooks and Monographs.
London: Chapman and Hall/CRC.
Lin, X., Lee, L.-f. (2010). GMM estimation of spatial autoregressive models with unknown heteroskedasticity. Journal of
Econometrics 157(1):34–52.
Liu, X., Lee, L.-f., Bollinger, C. R. (2010). An efficient GMM estimator of spatial autoregressive models. Journal of
Econometrics 159(2):303–319.
Newey, W. K., McFadden, D. (1994). Large Sample Estimation and Hypothesis Testing. In: Engle, R. F., McFadden, D. L., eds.
Handbook of Econometrics, Vol. 4. Elsevier, pp. 2111–2245.
Pace, R. K., Barry, R. (1997). Quick computation of spatial autoregressive estimators. Geographical Analysis 29(3):232–247.
954 S. TAŞPINAR ET AL.

Pinkse, J., Slade, M. E. (1998). Contracting in space: An application of spatial statistics to discrete-choice models. Journal
of Econometrics 85(1):125–154.
Prucha, I. R. (2014). Instrumental variables/method of moments estimation. In: Fischer, M. M., Nijkamp, P., eds. Handbook
of Regional Science. Berlin, Heidelberg: Springer, pp. 1597–1617.
Windmeijer, F. (2005). A finite sample correction for the variance of linear efficient two-step GMM estimators. Journal of
Econometrics 126(1):25–51.
Windmeijer, F. (2008). GMM for panel data count models. In: Mátyás, L., Sevestre, P., eds. The Econometrics of Panel Data.
Advanced Studies in Theoretical and Applied Econometrics, Vol. 46. Berlin, Heidelberg: Springer, pp. 603–624.