I need to estimate spatial econometric models with spatial lags of X (SLX) either alone, combined with spatial autoregressive models (SAR) or with spatial error model (SEM). When they are combined, they are called Spatial Durbin Model (SDM) or Spatial Durbin Error Model (SDEM), following Vega & Elhorst's (2015) paper "The SLX Model".
I intend to estimate all spatial panel models in R using splm package, which also requires spdep functions. In this sense, I created neighbors lists type Queen and k = 4 from a shape file:
> TCAL <- readOGR(dsn = ".", "Municipios_csv")
> coords <- coordinates(TCAL)
> contnbQueen <- poly2nb(TCAL, queen = TRUE)
> enter code herecontnbk4 <- knn2nb(knearneigh(coords, k = 4, RANN = FALSE))
Then I converted this neighbor list in a Weight matrix:
> W <- nb2listw(contnbk4, glist = NULL, style = "W")
attributes(W)
$names
[1] "style" "neighbours" "weights"
$class
[1] "listw" "nb"
$region.id
[1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"
[21] "21" "22" "23" "24" "25" "26" "27" "28" "29" "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "40"
[41] "41" "42" "43" "44" "45" "46" "47" "48" "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
[61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70" "71" "72" "73" "74" "75" "76" "77" "78" "79" "80"
[81] "81" "82" "83" "84" "85" "86" "87" "88" "89" "90" "91" "92" "93" "94" "95" "96" "97" "98" "99" "100"
[101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110" "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
[121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130" "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
[141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150" "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
[161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170" "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
[181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190" "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
[201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210" "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
[221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230" "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
[241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250" "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
[261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270" "271" "272" "273" "274" "275" "276"
$call
nb2listw(neighbours = contnbk7, glist = NULL, style = "W")
Next step, I created a formula for panel SAR and SEM models, which worked well and produced estimates:
> fmPanel <- Area ~ Dist + Land + CredAg
> vegSAR <- spml(fmPanel, data = veg, index = c("Mun","Year"), listw = W, model = "within", effect = "twoways", spatial.error = "none", lag = TRUE)
> vegSEM <- spml(fmPanel, data = veg, index = c("Mun","Year"), listw = W, model = "within", effect = "twoways", spatial.error = "b", lag = FALSE)
Then, I tried to estimate SLX, SDM and SDEM models by creating spatial lags of the covariates X:
> vegX <- pdata.frame(veg, index = c("Mun","Year")); class(vegX)
[1] "pdata.frame" "data.frame"
I then created pseries values:
> DistX <- vegX$Dist; class(DistX)
[1] "pseries" "numeric"
> LandX <- vegX$Land; class(LandX)
[1] "pseries" "numeric"
> CredAgX <- vegX$CredAg; class(CredAgX)
[1] "pseries" "numeric"
But the error happened when I applied slag function:
DistX <- slag(agSPX$Dist, listw = W)
Error in lag.listw(listw, xt) : object lengths differ
My panel data has 5 years and 276 regions. So, objects characteristics are:
> length(DistX)
[1] 1380
> length(W)
[1] 3
> length(W$weights)
[1] 276
So, I was wondering that, if I could transform W$weights in a matrix such as usaww used as example of slag function, I could apply function mat2listw and then use slag over X.
Could someone tell me where I am wrong?
You can also add the lagged explanatory variable to the data.
A reproducible example:
library(plm)
library(spatialreg)
library(splm)
# load data
data(Produc, package = "plm")
data(usaww, package = "splm")
d <- pdata.frame(Produc, index = c("state","year"), drop.index = FALSE)
# create a spatial explanatory variable
d$unemp_l <- slag(d$unemp, usaww)
# run model
m <- splm::spml(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp + unemp_l,
data = d, listw = mat2listw(usaww) , model="within")
summary(m)
Spatial panel fixed effects error model
Call:
splm::spml(formula = log(gsp) ~ log(pcap) + log(pc) + log(emp) +
unemp + unemp_l, data = d, listw = mat2listw(usaww), model = "within")
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-0.1211492 -0.0234013 -0.0040218 0.0167919 0.1787587
Spatial error parameter:
Estimate Std. Error t-value Pr(>|t|)
rho 0.542254 0.033772 16.056 < 2.2e-16 ***
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
log(pcap) 0.0090575 0.0251036 0.3608 0.71824
log(pc) 0.2152367 0.0234077 9.1951 < 2e-16 ***
log(emp) 0.7833003 0.0277672 28.2096 < 2e-16 ***
unemp -0.0014795 0.0011443 -1.2930 0.19603
unemp_l -0.0031210 0.0015790 -1.9766 0.04808 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
You can verify that the spatial explanatory is correct. For instance:
a <- usaww["ALABAMA",]
a <- a[a!=0]
a
FLORIDA GEORGIA MISSISSIPPI TENNESSE
0.25 0.25 0.25 0.25
mean(d[d$year=="1970" & d$state %in% names(a) , "unemp"])
[1] 4.525
d[d$state=="ALABAMA" & d$year=="1970", "unemp_l"]
ALABAMA-1970
4.525