Tutorial: Spatial econometrics with R 1

Spillover masking example : with a non spatial model

cbd <-  data.frame(y  =  c(42,  37,  30,  26,  30,  37,  42),
                   density = c(10, 20, 30, 50, 30, 20, 10),
                   distance = c(30, 20, 10, 0, 10, 20, 30),
                   row.names = paste("R",  1:7,  sep="")  ) 
head(cbd)

##     y density distance
## R1 42      10       30
## R2 37      20       20
## R3 30      30       10
## R4 26      50        0
## R5 30      30       10
## R6 37      20       20

Ordinary least squares model : \(Y=X \beta+\epsilon\) with \(\epsilon \sim \mathcal{N}\left(0, \sigma^{2} I_{n}\right)\)

obj.lm <- lm(y ~ density + distance-1, data = cbd)
summary(obj.lm)

## 
## Call:
## lm(formula = y ~ density + distance - 1, data = cbd)
## 
## Residuals:
##      R1      R2      R3      R4      R5      R6      R7 
## -0.9852  0.9926  0.9705 -1.5646  0.9705  0.9926 -0.9852 
## 
## Coefficients:
##          Estimate Std. Error t value Pr(>|t|)    
## density   0.55129    0.02064   26.71 1.37e-06 ***
## distance  1.24908    0.02839   43.99 1.15e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.284 on 5 degrees of freedom
## Multiple R-squared:  0.9991, Adjusted R-squared:  0.9987 
## F-statistic:  2649 on 2 and 5 DF,  p-value: 2.731e-08

Non-spatial prediction : \(\hat{Y}=X \hat{\beta}\) If we double the population density of region R2, only the travel time to the CBD for region R2 is affected:

cbd2=cbd 
cbd2$density[2]=2*cbd2$density[2]

yhat1 = predict(obj.lm, newdata = cbd) 
yhat2 = predict(obj.lm, newdata = cbd2) 
diff = yhat1-yhat2
diff 

##        R1        R2        R3        R4        R5        R6        R7 
##   0.00000 -11.02583   0.00000   0.00000   0.00000   0.00000   0.00000

Spillover masking example : with a spatial model

require("spdep")
require("spatialreg")

Consider the following spatial autoregressive model (LAG model) \(Y=\rho W y+X \beta+\epsilon\), with \(\epsilon \sim \mathcal{N}\left(0, \sigma^{2} I_{n}\right)\) We need to specify spatial neighbors through a weight matrix:

W.mat = matrix(c(0, 1, 0, 0, 0, 0, 0,
                 1, 0,  1,  0,  0,  0,  0,
                 0, 1,  0,  1,  0,  0,  0,
                 0, 0,  1,  0,  1,  0,  0,
                 0, 0,  0,  1,  0,  1,  0,
                 0, 0,  0,  0,  1,  0,  1,
                 0, 0,  0,  0,  0,  1,  0),7, 7,byrow=T)

W.listw <- mat2listw(W.mat)
Wnorm.listw = mat2listw(W.mat, style = "W")
Wnorm.mat = listw2mat(Wnorm.listw)

Z=rnorm(7,5,1) 
Zlag=Wnorm.mat%*%Z 
plot(Z,Zlag)

#changing the density in location 2 
cbd2=cbd 
cbd2$density[2]=2*cbd2$density[2]
yhat1.lm = predict(obj.lm, newdata = cbd) 
yhat2.lm = predict(obj.lm, newdata = cbd2) 
diff = yhat2 - yhat1
diff

##       R1       R2       R3       R4       R5       R6       R7 
##  0.00000 11.02583  0.00000  0.00000  0.00000  0.00000  0.00000

We can estimate the parameters :

obj.lag.norm <- lagsarlm(y ~ density + distance-1, 
                         data = cbd, listw=Wnorm.listw)
rho_hat = obj.lag.norm$rho

beta_hat = obj.lag.norm$coefficients
X = cbind(cbd$density,cbd$distance) 
X2 = cbind(cbd2$density,cbd2$distance)     
yhat1.lag = solve(diag(7) - rho_hat*Wnorm.mat)%*%X%*%beta_hat
yhat2.lag = solve(diag(7) - rho_hat*Wnorm.mat)%*%X2%*%beta_hat 
diff.lag = yhat2.lag-yhat1.lag
diff.lag  

##            [,1]
## [1,] 2.56991367
## [2,] 3.99688568
## [3,] 1.45570301
## [4,] 0.53110962
## [5,] 0.19632481
## [6,] 0.07956290
## [7,] 0.05115727

Spatial econometrics models allow to measure the effects of spillover

Illustration of estimation bias

library("raster")
library("cartography")
setwd("D:/Tap huan VIASM/Tinh hoa__VIASM_2021/Spatial econmetrics")
load("immob.Rdata")

immob.sp <- SpatialPoints(cbind(immob$longitude, immob$latitude))
immob.spdf <- SpatialPointsDataFrame(immob.sp, immob) 
proj4string(immob.spdf) <- CRS("+init=epsg:27572")

france_dep <- getData("GADM", country = "FRA", level = 2) 
france_dep <- spTransform(france_dep, CRS("+init=epsg:27572")) 
#save(france_dep, file = "france_dep.RData") 
#load("france_dep.RData")

#op <- par(oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0)) 
plot(france_dep, col = "#D1914D", border = "grey80") 
propSymbolsLayer(spdf = immob.spdf, var = "prix.vente", inches = 0.15)

#par(op)