40 Spatial econometrics introduction
40.1 Why do we need spatial econometrics?
‘Everything is related to everything else, but near things are more related than distant things’ (Tobler 1970). This is the first law of geography. This is the reason why we need spatial econometrics. The spatial econometrics is to account for the spatial dependence in the data. The spatial dependence is the correlation between the observations that are close to each other. This is a common feature in many datasets. For example, in the housing market, the price of a house is related to the prices of the houses nearby. In the crime data, the crime rate in one neighborhood is related to the crime rate in the neighboring neighborhoods. In the health data, the health outcomes of people in one area are related to the health outcomes of people in the neighboring areas. In the environmental data, the pollution level in one area is related to the pollution level in the neighboring areas. In the economic data, the economic performance of one region is related to the economic performance of the neighboring regions. In the social network data, the behavior of one person is related to the behavior of the people in the neighboring areas.
Econometricians love OLS. It is worth being precise about what spatial dependence breaks. Spatially correlated errors alone (the term \(u = \rho W u + \epsilon\) below, with exogenous regressors) do not make the OLS coefficients inconsistent — they make the conventional standard errors wrong, so inference is invalid even though the point estimates are fine. OLS coefficients become biased/inconsistent in the genuinely simultaneous or omitted-variable cases: a spatially lagged outcome \(Wy\) on the right-hand side (endogenous by construction), omitted spatially correlated confounders, or endogenous spillovers. The sections below treat these cases separately: S2SLS for the spatial-lag-of-outcome model, and spatial-error models for the correlated-error case.
40.2 A general linear spatial model
\[ y = \lambda W y + X \beta + W X \gamma + u \]
\[ u = \rho W u + \epsilon \]
where \(Y\) is the dependent variable, \(X\) is the matrix of independent variables, \(\beta\) is the vector of coefficients, and \(u\) is the error term. The error term \(u\) is spatially correlated.
We can see the major difference is \(W\). The term \(W y\) is called “spatially lagged y”, which in my opinion is misleading. This term is likely from time-series counter-part of lagged \(y\), but there is nothing about lags in the spatial data. \(W\) is the spatial weights matrix. It is a matrix that describes the spatial structure of the data. The simplest example of \(W\) is the case that we define an element of \(W\), \(w_{i,j}\) to be 1 if unit \(i\) and \(j\) are neighbors, 0 if not.
In this simplest case, our model is saying we have a pattern we need to take account of: neighboring units are related, they seem to display similar pattern in terms of our study. This commonality of neighbors can be causal, for example, if we study the effect of minimum wage effect on employment, then if the neighboring state’s wage increase might negatively impact the focal state’s employment, because people go for higher wages. This is the spatial spillover effect. This can also be non-causal, for example, if we study the effect of temperature on crop yield, then the neighboring state’s temperature might be a good predictor of the focal state’s temperature. This is the spatial autocorrelation. In either case, we need to take account of the spatial structure in the data.
The reason behind the spatial pattern can be all kinds. It can be \(x_i\) causing \(y_j\), or \(y_i\) causing \(y_j\), or some unobserved factors causing both \(y_i\) and \(y_j\), or some other reasons. Ignoring this spatial structure could make OLS biased. Betz et al (2020) show that ignoring spatial interdependence in the outcome results in asymptotically biased estimates even when instruments are randomly assigned.
40.3 Estimation with S2SLS for a less general model
Because we have the term \(W y\), we introduced endogeneity. \(W y\) is necessarily correlated with \(u\), even if \(u\) is not spatially correlated. Suppose \(u\) is not spatially correlated.
The spatial two stage least squares (S2SLS) estimator is one solution. Usually we use \(X\) and spatial lags of \(X\) as instruments. Typically instruments \(H\) is \((X, W X, W^2 X)\).
S2SLS is similar to 2SLS. The first stage is \[ \hat Z = PZ = [X, WX, \hat {W y}] \] where \[\hat {W y} = PWy\] and \[P = H(H'H)^{-1}H'\].
Then \[ \hat \theta_{S2SLS} = (\hat Z'Z)^{-1} \hat Z'y \] where \(\hat \theta_{S2SLS}\) is the full coefficient vector \((\beta, \gamma, \lambda)\).
The first stage is to regress \(WY\) on \(H\) . The second stage is to regress \(Y\) on the predicted \(WY\) from the first stage and \(X\) and \(W X\). The coefficient of \(X\) is the S2SLS estimator.
40.4 Estimation with a more general model
Suppose instead the error term is spatially correlated
\[ u = \rho W u + \epsilon \]
Then the estimation is more complicated. The generalized spatial two stage least squares (GS2SLS) estimator is needed. See Kelejian and Prucha (1998) for details.
40.5 an example
This example is from R package “sphet”.
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 avail- able, 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 accessibil- ity 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 spatial weights matrix is a sphere of influence neighbors list also available from spdep once the Boston data are loaded. Here the weight matrix comes with the data. Usually it’s the most heavy lifting part of the study.
The function read.gwt2dist reads a ‘.GWT’ file (e.g., generated using the function distance). In this example we are using the matrix of tract point coordinates projected to UTM zone 19 boston.utm available from spdep to generate a ‘.GWT’ file of the 10 nearest neighbors.
id1 <- seq(1, nrow(boston.utm))
tmp <- distance(boston.utm, region.id = id1, output = TRUE,
type = "NN", nn = 10, shape.name = "shapefile",
region.id.name = "id1", firstline = TRUE,
file.name = "boston_nn_10.GWT")
coldist <- read.gwt2dist(file = "boston_nn_10.GWT", region.id = id1, skip = 1)
class(coldist)[1] "sphet" "distance" "nb" "GWT" summary(coldist)
Number of observations:
n: 506
Distance summary:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.5441 0.9588 1.5843 2.0848 2.6389 11.6388
Neighbors summary:
Min. 1st Qu. Median Mean 3rd Qu. Max.
10 10 10 10 10 10 Spatial models in sphet are fitted by the functions stslshac and gstslshet. “stslshac” is Spatial two stages least square with HAC standard errors. “gstslshet” is GM estimation of a Cliff-Ord type model with Heteroskedastic Innovations, which is the general model we discussed before, with spatially correlated errors.
The theoretical model above is the Spatial Durbin Model (SDM), which includes \(WX\) as regressors (with their own coefficient \(\gamma\)), not just as instruments for \(Wy\). If we only pass \(X\) (no \(WX\)) in the formula below, stslshac’s W2X argument still uses spatial lags of \(X\) internally, but only as instruments for the endogenous \(Wy\) – \(\gamma\) is never estimated, so the fitted model is actually the simpler Spatial Autoregressive (SAR) model (\(\gamma=0\)), not SDM. To match the SDM specification, include \(WX\) explicitly as regressors:
boston_x_vars <- c("CRIM","ZN","INDUS","CHAS","NOX","RM","AGE","DIS","RAD","TAX","PTRATIO","B","LSTAT")
boston.c2 <- boston.c
for (v in boston_x_vars) {
xv <- boston.c[[v]]
if (is.factor(xv)) xv <- as.numeric(as.character(xv))
boston.c2[[paste0("W_", v)]] <- lag.listw(listw, xv)
}
sdm_formula <- as.formula(paste(
"log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT) +",
paste0("W_", boston_x_vars, collapse = " + ")
))
res <- stslshac(sdm_formula,
data = boston.c2, listw,
distance = coldist, type = "Triangular", HAC = TRUE)
summary(res)======================================================
======================================================
Spatial Lag Model
HAC standard errors
======================================================
======================================================
Call:
stslshac(formula = sdm_formula, data = boston.c2, listw = listw,
HAC = TRUE, distance = coldist, type = "Triangular")
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-0.71074 -0.06258 -0.00291 0.06279 0.66519
Coefficients:
Estimate SHAC St.Er. t-value Pr(>|t|)
Wy 0.71491435 0.09797746 7.2967 2.949e-13 ***
(Intercept) 1.74597474 0.46867690 3.7253 0.0001951 ***
CRIM -0.00548837 0.00153211 -3.5822 0.0003407 ***
ZN 0.00081016 0.00045400 1.7845 0.0743454 .
INDUS -0.00070720 0.00285254 -0.2479 0.8041961
CHAS1 -0.05854830 0.03984589 -1.4694 0.1417328
I(NOX^2) -0.03240666 0.21228222 -0.1527 0.8786677
I(RM^2) 0.00846160 0.00289595 2.9219 0.0034793 **
AGE -0.00150525 0.00047265 -3.1847 0.0014489 **
log(DIS) -0.16117615 0.07522729 -2.1425 0.0321515 *
log(RAD) 0.04219650 0.01648649 2.5595 0.0104835 *
TAX -0.00045089 0.00012647 -3.5652 0.0003636 ***
PTRATIO -0.01336079 0.00581140 -2.2991 0.0215013 *
B 0.00058124 0.00011965 4.8578 1.187e-06 ***
log(LSTAT) -0.21713640 0.03831649 -5.6669 1.454e-08 ***
W_CRIM -0.00451534 0.00300618 -1.5020 0.1330927
W_ZN -0.00133281 0.00058617 -2.2738 0.0229807 *
W_INDUS 0.00079183 0.00351396 0.2253 0.8217152
W_CHAS 0.09604488 0.06127241 1.5675 0.1169964
W_NOX -0.34866680 0.32382701 -1.0767 0.2816112
W_RM -0.06511363 0.05723288 -1.1377 0.2552474
W_AGE 0.00175310 0.00062322 2.8130 0.0049086 **
W_DIS 0.01948831 0.01383400 1.4087 0.1589164
W_RAD 0.00203459 0.00345353 0.5891 0.5557722
W_TAX 0.00028960 0.00019481 1.4866 0.1371284
W_PTRATIO 0.00732717 0.00807107 0.9078 0.3639669
W_B -0.00054543 0.00015750 -3.4630 0.0005342 ***
W_LSTAT 0.00895071 0.00427009 2.0961 0.0360698 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The function gstslshet is used to estimate the GS2SLS model. Unlike stslshac, gstslshet’s GMM machinery hits an exactly-singular matrix on this dataset once the \(WX\) regressors are added (likely due to near-collinearity introduced by the additional lagged terms combined with its own internal instrument construction), so we keep the original (SAR, no \(WX\)) specification here — this estimates the spatial-error model with \(\gamma=0\), not the full Durbin extension shown above for stslshac.
Call:
gstslshet(formula = 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 = listw, initial.value = 0.2)
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.569385 -0.073165 -0.001675 0.000531 0.071500 0.740309
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 2.51316605 0.26749367 9.3952 < 2.2e-16 ***
CRIM -0.00662744 0.00144522 -4.5858 4.523e-06 ***
ZN 0.00038299 0.00036563 1.0475 0.2948733
INDUS 0.00159352 0.00179772 0.8864 0.3753967
CHAS1 -0.00447974 0.03689065 -0.1214 0.9033480
I(NOX^2) -0.27295899 0.11561412 -2.3609 0.0182283 *
I(RM^2) 0.00744059 0.00199637 3.7270 0.0001937 ***
AGE -0.00045400 0.00045572 -0.9962 0.3191333
log(DIS) -0.16517174 0.03484858 -4.7397 2.140e-06 ***
log(RAD) 0.07453521 0.01752830 4.2523 2.116e-05 ***
TAX -0.00041956 0.00010763 -3.8981 9.697e-05 ***
PTRATIO -0.01412661 0.00410143 -3.4443 0.0005725 ***
B 0.00035970 0.00011182 3.2168 0.0012964 **
log(LSTAT) -0.24593826 0.03213364 -7.6536 1.954e-14 ***
lambda 0.42407826 0.04463747 9.5005 < 2.2e-16 ***
rho 0.29587455 0.08614291 3.4347 0.0005932 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In stata, you would use “spivregress” which replaced user-written “spivreg”. Stata has a set of commands for spatial econometrics. Check out stata’s “sp” manual.