sphet: Spatial Models with Heteroskedastic Innovations in R

sphet is a package for estimating and testing spatial models with heteroskedastic innovations. We implement recent generalized moments estimators and semiparametric methods for the estimation of the coeﬃcients variance-covariance matrix. This paper is a general description of sphet and all functionalities are illustrated by application to the popular Boston housing dataset. The package in its current version is limited to the estimators based on Arraiz, Drukker, Kelejian, and Prucha (2010); Kelejian and Prucha (2007, 2010). The estimation functions implemented in sphet are able to deal with virtually any sample size.


Introduction
sphet is a package for estimating and testing a variety of spatial models with heteroskedastic innovations.The estimation procedures are based on generalized moments (GM).
An increasing number of datasets contain information on the location of the observations along with the values of the variables of interest.Taking into account the spatial nature of the data can improve efficiency or, in some cases, even be essential for consistency.This increasing interest in spatial data is corroborated by the large number of packages in the R language (R Development Core Team 2010) for analyzing multiple typologies of spatial data under different methodological perspectives.Among these alternatives, spdep is one of several packages that deals with spatial dependence (Bivand 2001(Bivand , 2002(Bivand , 2006;;Bivand and Gebhardt 2000;Bivand and Portnov 2004).spdep includes functions to create and manipulate spatial objects (i.e., creating spatial weights matrices).It also contains a collection of tests for spatial autocorrelation and functions for estimating spatial models.For what concerns estimation features, general method of moments (GMM), instrumental variables (IV) and maximum likelihood (ML) are supported.sphet complements but does not overlap with the econometric features already available in spdep.Specifically, spdep focuses on spatial lag and error models, whereas our departure point is the following model: u = ρW u + ε Kelejian and Prucha (2007) do not assume any specific structure for the disturbance process.
Their focus is to develop an estimation theory (for the regression parameters) that is robust against possible misspecification of the disturbances.Nonetheless, the general assumptions made on the disturbance process cover the spatial autoregressive process as a special case.
Although Model (1) is well established in the econometrics literature, regional scientists and geographers generally follow a different approach.They start from an OLS regression and try to determine (via Lagrange multipliers (LM) tests on estimated residuals) whether the true data generating process is a spatial error, or a spatial lag model.This is unfortunate because the spatial patterns implied by Model (1) are richer than those implied by either the spatial error or the spatial lag model (Kelejian and Prucha 1998).We believe that providing this alternative implementation is a useful contribution to both scientific communities.Related to this, Kelejian and Prucha (2010) give results concerning the joint asymptotic distribution of IV and GMM estimators for Model (1).Their results permit testing the (joint) hypothesis of no spatial spillovers originated from the endogenous variables or from the disturbances.As a consequence, if one of the corresponding coefficients turns out not to be statistically different from zero, one could still go back to the estimation of a reduced model.
In sphet, we only concentrate on GM and IV methods, leaving aside the ML approach.
Reasons for this will be discussed throughout the paper.In general terms, GM requires weaker assumptions than ML.Additionally, there are still various unsolved problems related to the ML approach.Numerical difficulties related to the computation of the Jacobian term (Kelejian andPrucha 1998, 1999) may potentially limit the application of ML to large datasets.However, various solutions have been proposed in the literature that have attenuated the problem (see e.g., Ord 1975;Smirnov and Anselin 2001;Pace 1997;Pace and Barry 1997;Barry and Pace 1999;Pace and LeSage 2004, among others).Lee (2004) derives the conditions that ensure consistency and asymptotic normality of ML estimators for the general spatial model considered but some of the assumptions are stronger than those required by GM.On the other hand, there is reasonably general theory for the GM approach in the cross sectional case (Kelejian and Prucha 1998, 1999, 2010).
One final point relates to the possible presence of heteroskedasticity in the innovations of the model.Since spatial units may differ in important characteristics (e.g., size) homoskedasticity is a strong assumption that may not hold in many applied spatial problems.Anselin and Lozano-Gracia (2008) and Baltagi, Egger, and Pfaffermayr (2008) are two typical examples of empirical applications that require the use of spatial heteroskedasticity and autocorrelation consistent (HAC) estimators.The estimation theory developed by Lee (2004) for the quasimaximum likelihood estimator under the assumption of homoskedastic innovations does not carry over to the case of heteroskedastic innovations.To support this, Arraiz et al. (2010) provide simulation evidence that when the innovations are heteroskedastic, ML produces inconsistent estimates.We implement various GM and IV procedures to obtain spatial HAC estimators of the variance-covariance matrix of the model coefficients.In its current version, the package is limited to methodologies implemented in Kelejian andPrucha (2007, 2010) and Arraiz et al. (2010).The gstslshet code was tested against the original Stata (StataCorp 2007) code used to produce the simulation results in Arraiz et al. (2010).The function stslshac give results like those implemented in Python (van Rossum 1995) and used in Anselin and Lozano-Gracia (2008).Unfortunately, none of the cited implementations are available to the general public.The estimation functions implemented in sphet are able to handle virtually any sample size.
To overcome part of the technical difficulties that arise in the implementation, we make extensive use of classes of objects defined in spdep as for example spatial weights matrix objects (of class listw).Furthermore, we also make substantial use of code from the Matrix package (Bates and Maechler 2010).
The remainder of the present paper is a general description of sphet and all functionalities are illustrated by application to the popular Boston housing dataset (Harrison and Rubinfeld 1978).
The Boston data contain information on property values and characteristics in the area of Boston, Massachusetts and has been widely used for illustrating spatial models.Specifically, there is a total of 506 units of observation for each of which a variety of attributes are available, such as: (corrected) median values of owner-occupied homes (CMEDV); per-capita crime rate (CRIM); nitric oxides concentration (NOX); average number of rooms (RM); proportions of residential land zoned for lots over 25,000 sq.ft (ZN); proportions of non-retail business acres per town (INDUS); Charles River dummy variable (CHAS); proportions of units built prior to 1940 (AGE); weighted distances to five Boston employment centers (DIS); index of accessibility to highways (RAD); property-tax rate (TAX); pupil-teacher ratios (PTRATIO); proportion of blacks (B); and % of the lower status of the population (LSTAT).
The dataset with Pace's tract coordinates is available to the R community as part of spdep.
The sphet package is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=sphet.After installation, it can be loaded via library("sphet") and data("boston", package = "spdep") loads the Boston data from spdep.

Tools
sphet supports a series of capabilities to generate distance objects that are required by the semiparametric estimation of the variance-covariance matrix discussed in Section 3.These functionalities can be accessed through the functions distance and read.gwt2dist.
The function distance reads points coordinates and generates a matrix.The object created is similar to the content of a '.GWT' file.A '.GWT' file format is a well known structure among spatial statisticians and econometricians.These (memory efficient) types of weights matrices can be calculated using GeoDa (Anselin, Syabri, and Kho 2006) and other software.
In fact, one advantage of R over other packages is the availability of these models that are well established among users.The generated object is made up of three columns.The third column lists the distances, while the first two columns contain the id of the two points for which the distance is calculated.
Currently, five distance measures are implemented, namely: and the Great Circle distance.For details on how to calculate this last distance measure one should see the function rdist.earth in the package fields (Furrer, Nychka, and Sain 2009).
The following instructions demonstrate the usage of the function distance.We first generate a set of XY-coordinates corresponding to one hundred points.The first coordinate is a random sample from the uniform distribution on the interval (0, 70), whereas the second coordinate is generated on the interval (−30, 20).The object coord1 is a matrix whose first column is intended to contain the identification for the points.
R> set.seed(1234)R> X <-runif(100, 0, 70) R> Y <-runif(100, -30, 20) R> coord1 <-cbind(seq(1, 100), X, Y) The (optional) id variable has the principal scope of providing the ordering of the observations.When specified, it could be the first column of the argument coord.Alternatively, it could be specified separately as region.id.When region.id is not NULL and coord has three columns (i.e., an id variable has been specified twice) the function performs some checks to make sure that the two variables point to the same ordering.On the other hand, if an id variable is not specified at all, it is assumed to be a sequence from one to the number of observations (i.e., the number of coordinates).R> thm1 <-distance(coord = coord1, region.id= NULL, output = FALSE, + type = "inverse", measure = "euclidean") R> print(thm1[1:15, ]) The measure argument specifies the distance measure and should be one of "euclidean", "gcircle", "chebyshev", "braycur", and "canberra".The type argument is used to define the distance criteria and should be one of "inverse", "NN" or "distance".Both "inverse" and "distance" can be specified along with a cutoff.The cutoff takes up three values: 1, 2, and 3 indicating the lower, median and upper quantile of the distance distribution.Specifically, when cutoff is set to 1, only observations within a distance less than the first quantile are neighbors to each other.All other interactions are considered negligible."NN" (nearest neighbors) should be specified along with nn, the argument to define the number of nearest neighbors, as it is illustrated in the following example.When output is TRUE, the function writes the data to a file.The output file can have any format.In particular, it could be a '.GWT' file.When firstline is TRUE, an header line is added to the '.GWT' file.The first element is simply a place holder, the second is the number of observations.The name of the shape file and of the id variable can be specified by the options shape.nameand region.id.namerespectively.If an output file is produced, the name of the file can be set by file.name.
In this example we are using the matrix of tract point coordinates projected to UTM zone 19 boston.utmavailable from spdep to generate a '.GWT' file of the 10 nearest neighbors.The file 'boston_nn_10.GWT" is then inputted to the function read.gwt2distalong with the region.id.
It is worth noticing that the function read.gwt2distcould also read other extensions (such as '.txt').It is important, however, that the input file exhibits the general format described above.When the file has a '.GWT' extension, the number of observations is generally retrieved from the first line.Alternatively, it is fixed to the length of the (unique) region.idvariable.
The argument skip determines the number of lines to disregard before reading the data.The value is an object of class distance.We generate a new class of objects to be able to perform some of the checks necessary to make sure that the distance measure specified in the function stslshac is appropriate.

R> class(coldist)
[1] "sphet" "distance" "nb" "GWT" A summary method is available to print some of the basic features of the distance object.
In particular, the total number of observations and some general descriptive statistics for both distances and neighbors are displayed.We believe that this information is of guidance while choosing the type of bandwidth to employ in the spatial HAC estimation discussed in Section 3.

Estimation functions
Spatial models in sphet are fitted by the functions stslshac and gstslshet.Below, we first review some of the theory and then demonstrate the use of these functions.

Spatial two stages least squares with HAC standard errors
Consider the following spatial model: or also, more compactly . The presence of the spatially lagged dependent variable W y introduces a form of endogeneity.Under typical specifications, W y will be correlated with the disturbances ε, which motivates an instrumental variable approach.The spatial two stage least squares (S2SLS) regression is a straightforward extension of the "classical" two stage least squares procedure.The selection of instruments as an approximation to ideal instruments has been extensively discussed in the literature (see e.g., Kelejian andPrucha 1998, 1999;Kelejian, Prucha, and Yuzefovich 2004;Lee 2003, among others) and is based on the reduced form of the model.In empirical problems, the matrix of instruments can be defined in the following way: where, typically, q ≤ 2. The matrix of instruments implemented in sphet is The S2SLS estimator for the parameter vector γ can be obtained as: where Ẑ = P Z = [X, W y], W y = P W y and P = H(H H) −1 H . Statistical inference is generally based on the asymptotic variance covariance matrix: with σ2 = e e/n and e = y − Z γS2SLS .Kelejian and Prucha (2007) propose a HAC estimation of the variance covariance matrix of Model (2).The spatial HAC estimator is robust against possible misspecification of the disturbances and allows for (unknown) forms of heteroskedasticity and correlation across spatial units.The disturbance vector is assumed to be generated by the following very general process: where ξ is a vector of innovations and R is an n × n non stochastic matrix whose elements are not known.Additionally, R is non-singular and the row and column sums of R and R −1 are bounded uniformly in absolute value by some constant (for technical details see Lee 2002Lee , 2003Lee , 2004;;Kelejian and Prucha 1998, 1999, 2004).This specification of the error term covers SARMA(p,q) processes as special cases.Even if we assume such a general specification for the disturbance process we still have to be concerned about possible misspecifications (e.g., due to an incorrect specification of the weights matrices).The asymptotic distribution of corresponding IV estimators involves the variance covariance matrix where Σ = RR denotes the variance covariance matrix of ξ.Kelejian and Prucha (2007), propose to estimate the r, s elements of Ψ by: where the subscripts refer to the elements of the matrix of instruments H and the vector of estimated residuals ε.Kelejian and Prucha (2007) also contains a generalization to several distance measures.In that case, the expression for the spatial HAC estimator of the true variance covariance matrix assumes a slightly different form.However, we only implement the estimator based on the case of a single measure.K() is a Kernel function used to form the weights for the different covariance elements.The Kernel function is a real, continuous and symmetric function that determines the pairs of observations included in the cross products in (9).The Kernel function is defined in terms of a distance measure.More specifically, d * ij represent the distance between observations i and j and d is the bandwidth.Note that Kelejian and Prucha (2007) allows for the case where the researcher measures these distances with error.More in detail, the distance measure employed by the researcher is given by where υ ij denotes the measurement error.The only assumption made on the random measurement error is that it is independent on the innovations of the model ξ.
In other words, the bandwidth plays the same role as in the time series literature; Together with the Kernel function it limits the number of sample covariances.Furthermore, the bandwidth can be assumed either fixed or variable.A fixed bandwidth corresponds to a constant distance for all spatial units.On the other hand, a variable bandwidth varies for each observation (i.e., the distance corresponding to the n-nearest neighbors).
Based on the spatial HAC estimator of Ψ given in (9), the asymptotic variance covariance matrix ( Φ) of the S2SLS estimator of the parameters vector is given by: Therefore, small sample inference can be based on the approximation γ ∼ N (γ, n −1 Φ).

Demonstration
The function that deals with the spatial HAC estimator is stslshac.Crucial arguments are listw, distance, type and bandwidth.stslshac requires the specification of two different lists: one for the spatial weights matrix W and one to define the distance measure d.As in spdep, listw is the argument that handles the spatial weights matrix W in the form of a list.The object listw can be generated for example by the function nb2listw available in spdep.On the other hand, the argument distance that specifies the distance measure, is an object of class distance created for example by read.gwt2dist.Note that the two objects, although belonging to a different class, may be generated according to the same definition.The argument type deals with the specification of the kernel function.Currently, six different kernels are available: − cos(6πz)/5)).
If the kernel type is not one of the six implemented, the function will terminate with an error message.Note that this last result corresponds to the one obtained by stsls in spdep (with robust = FALSE).Both functions have a logical argument (W2X) that if set to TRUE (the default) uses the matrix of instruments H = (X, W X, W 2 X) in the spatial two stages least squares.Since (Kelejian et al. 2004) show the advantages of including W 2 X in the matrix of instruments, we strongly recommend to leave the argument W2X at its default value.In this example, no substantial differences are observed in terms of significance of the coefficients when using the robust estimator.It would be good practice to always estimate HAC standard errors at least to compare them with traditional results.If this leads to different significance levels, one should always present robust results.Residual variance (sigma squared): 0.020054, (sigma: 0.14161) The heteroskedasticity correction leads to results similar to those obtained with the spatial HAC estimator implemented in sphet (on the Boston housing data).The spatial HAC seems to be more conservative in that it changes the significance level of some of the variables (i.e., NOX2, RM2 and B).
The argument bandwidth by default sets the bandwidth for each observation to the maximum distance for that observation.Alternatively, a fixed bandwidth can be used as in the next example that fixes as bandwidth the maximum distance (overall).0.00028845 0.00016362 1.7629 0.0779146 .log(LSTAT) -0.23984212 0.03955045 -6.0642 1.326e-09 *** ---Signif.codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 summary and print methods are available.The summary method provides a short description of the model followed by a printout of the most recent call, a summary of the residuals, and the table of the results.The first row of the table of estimated coefficients is always the spatial lag of the dependent variable W y. Note that the name of the column of the standard errors clearly make reference to the use of the spatial HAC estimator.On the other hand, the print method simply prints basic information.
One final point deserves to be mentioned.As we anticipated in the introduction, we are not giving much attention to the ML approach because of various shortcomings.It is technically possible, however, to make heteroskedasticity corrections to standard errors in a ML context using functions in the lmtest (Zeileis and Hothorn 2002) and sandwich (Zeileis 2004(Zeileis , 2006) ) packages.1An example on how to perform the correction is given in the help file of the function bptest.sarlmavailable from spdep.We run this example on the Boston data set and we did not observe big shifts in the estimated coefficients or substantial differences in their significance.

General spatial two stage least squares
In Section 3.1 we assumed a very general form for the disturbance process.As an alternative, one could assume that the disturbance process is known to follow a first order spatial autoregressive process where the innovations ξ 1 , . . ., ξ n are assumed independent with zero mean and non-constant variance σ 2 i .Kelejian and Prucha (2010) define a GM estimator for ρ and give results for both consistency and asymptotic normality.Note that Kelejian and Prucha (1999) had only proven consistency of the GM estimator under the assumption of homoskedastic innovations.In other words, ρ was seen as a nuisance parameter (e.g., the distribution of the regression's parameters does not depend on the estimator for ρ).Kelejian and Prucha (2010) also give results for the joint asymptotic distribution of the GM estimator for ρ and the IV estimators for γ.The suggested estimation procedure consists of two steps alternating GM and IV estimators.Each of the two steps includes sub-steps.In the first step, γ is estimated by S2SLS with the matrix of instruments defined in Equation 4. The estimated residuals from the first step are employed to obtain a sample analogue of the moment conditions (for greater detail on the specification of the moments conditions see Kelejian and Prucha 2010;Arraiz et al. 2010).An initial GM estimator for ρ is defined in terms of these moment conditions and S2SLS residuals.The initial estimator can be viewed as an unweighted nonlinear least square estimator.Although fully consistent, it is not efficient because of a lack of weighting.This is why in the third sub-step of the first step, an efficient GM estimator for ρ is derived based on S2SLS residuals and moments conditions appropriately weighted by an estimator of the variance covariance matrix of the limiting distribution of the normalized sample moments.Following Kelejian and Prucha (1998), the consistent and efficient estimator for ρ is used to take a spatial Cochrane-Orcutt transformation of the model.In the second step of the estimation procedure the transformed model is estimated by S2SLS: this is the generalized spatial two-stage least square (GS2SLS) procedure presented in Kelejian and Prucha (1998).Specifically, the GS2SLS estimator for the parameter vector γ is defined as: where y * = y − ρW y, Z * = Z − ρW Z, Ẑ * = P Z * , and P = H(H H) −1 H .In the second and final sub-step of the second step, new sample moments are obtained by replacing S2SLS residuals by GS2SLS residuals obtained from Equation 12.The efficient GM estimator for ρ based on GS2SLS residuals is obtained from where the weighting matrix Υ−1 is an estimator of the variance covariance matrix of the limiting distribution of the sample moments.Under the assumptions made in Kelejian and Prucha (2010), γ and ρ are both consistent and asymptotically (jointly) normal.Therefore, small sample inference can be based on the following estimator for the variance covariance matrix: Finally, P , J, Σ and a are all expressions (based on the instruments and the transformed variables) needed to estimate the asymptotic variance covariance matrix of the moment conditions.The tilde explicitly denotes that all quantities are based on estimated residuals from GS2SLS procedure.
As a final point, a joint test of significance on the spatial coefficients can be derived.Let B be a 1 × (K + 2) matrix and θ = [γ, ρ ] the (K + 2) × 1 vector of estimated parameters (including the two spatial parameters).A restriction (e.g., that both spatial parameters are zero) can then be formulated in terms of Bθ = 0 A Wald test can then be based on (Greene 2003): and under H 0 will have a chi-squared distribution with one degree of freedom (i.e., the number of rows in B).
A complication appears from a computational perspective.The estimation of the variance covariance matrix of the limiting distribution of the normalized sample moments based on two stages least squares residuals involves the inversion of an n × n matrix.This is the main reason for transforming the object of class listw into a sparse Matrix and use code from the Matrix package to calculate the inverse.However, for very large problems the inversion could still be computationally intensive.In these cases the inverse can be calculated using the approximation where the last element of the sum depends on each specific W and ρ and is determined through a condition. 2For particular spatial weights matrices the results obtained using the approximation could be very close to the actual inverse of the variance covariance matrix.Furthermore, the inverse only influences the expression of the estimated variance covariance matrix of the limiting distribution of the normalized sample moments based on 2SLS residuals.In other words, small differences in the weighting matrix may imply even smaller differences in the estimated value of the spatial parameter resulting from the optimization procedure.As an example, on the Boston data the value of ρ resulting from the correct inverse is 0.1778462.If using the the approximation the value of ρ turns out to be 0.1779276.We would suspect that with larger datasets (for sparse W ) the difference should be even smaller.
A second issue is related to the initial values of the optimization search for the parameter ρ. 3 The default is to start from 0.2.As an alternative the user can either specify a different value or take as initial value the estimated coefficient of a regression of the S2SLS residuals on their spatial lag (i.e., fit a spatial autoregressive model on the residuals).The initial value in successive steps is the optimal parameter in previous steps.

Demonstration
The function that allows estimating the model described in this Section is gstslshet.It is also possible to estimate a restricted version of the model for which the parameter λ is set to zero by changing the logical argument sarar.Such a model is generally referred to in the literature as a spatial error model (Anselin 1988).The syntax of the function is straightforward in its basic arguments.The model to be estimated is described by a formula, an optional data.framecan be specified, and the spatial weights matrix is an object of class listw.The argument initial.valuemanages the starting point of the optimization process in the search for the optimal ρ.The default value for initial.value is 0.2.Any other numeric value (within the search interval) is acceptable.Alternatively, if initial.value is set to "SAR" the optimization will start from the estimated coefficient of a regression of the 2SLS residuals over their spatial lag.The argument inverse (default: TRUE) should be altered only when strictly necessary depending on the number of cross sectional observations.In this case, the inverse will be calculated by Equation 15.The precision of the approximation can be managed through the argument eps.The next example illustrates the use of the approximated inverse in the context of a model where λ is assumed to be zero (sarar = FALSE).Max.-0.832000 -0.090200 -0.002940 0.000891 0.095800 0.890000 nium Science Initiative Program "Regional Sciences and Public Policy" (MIDEPLAN -Chile).Usual disclaimers apply.
R> res <-stslshac(log(CMEDV) ~CRIM + ZN + INDUS + CHAS + I(NOX^2) + + I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT), + data = boston.c,listw, distance = coldist, HAC = FALSE) R> summary(res) stsls allows an heteroskedasticity correction to the coefficients' variance covariance matrix by setting the argument robust to TRUE.The additional argument legacy chooses between two different implementations of the robustness correction.When it is set to FALSE (the default used in our examples), a White consistent estimator of the variance-covariance matrix is provided.On the other hand, if legacy equals TRUE a GLS estimator is performed that yields different coefficient estimates.Results are displayed in the following example.