library(SDPDmod)
This vignette gives a few examples on how to use the
blmpSDPD
and SDPDmod
functions from the
SDPDmod
package.
The general spatial static panel model takes the form:
\[\begin{equation} \begin{aligned} y_{t}&=\rho W y_{t} + X_{t} \beta + W X_{t} \theta + u_{t}, \\ u_{t} &=\lambda W u_{t}+\epsilon_{t} \end{aligned} \label{eq:mod1}\tag{1} \end{equation}\]where the \(N \times 1\) vector \(y_{t}\) is the dependent variable, \(X_{t}\) is a \(N \times k\) matrix of \(k\) explanatory variables and \(W\) is a spatial weights matrix. \(N\) represents the number of units and \(t=1,...,T\) are the time points. The spatial lags for the vector of covariates is denoted with \(WX_t\). Spatial interaction in the error term \(u_{it}\) is included with the \(\lambda\) coefficient and it is assumed that \(\epsilon_t\) is independently and identically distributed error term for all \(t\) with zero mean and variance \(\sigma^2\). \(\rho\) is the spatial autoregressive coefficient, \(\lambda\) the spatial autocorrelation coefficient, \(\beta\) is a vector of response parameters for the covariates.
Figure 1. The different spatial panel models and their dependencies.
Note: SAC - spatial autoregressive combined,
SDM - spatial Durbin model, SDEM - spatial Durbin error model, SAR -
spatial autogregressive model (or Spatial lag model), SEM - spatial
error model, SLX - spatially lagged X model.
However, model (1) suffers from identification problems (Elhorst 2010). Figure 1 shows the identifiable models, which are derived by imposing some restrictions on model (1). If all spatial coefficients are zero (\(\rho = \theta = \lambda = 0\)), then the corresponding model will be the ordinary least-squares model (OLS). If in each of the models in figure (1) \(\eta W y_{t-1} +\tau y_{t-1}\) is added, we get the corresponding dynamic panel data models. \(y_{(t-1)}\) denotes the time lag of the dependent variable and \(Wy_{(t-1)}\) is the time-space lag. \(\tau\) and \(\eta\) are the response parameters for the lagged variable.
LeSage and Parent (2007) and LeSage (2014) suggest a Bayesian approach for
selecting an appropriate model. The calculation of the posterior
probabilities for 6 models (SAR, SEM, SDM, SDEM, SLX, OLS) is possible
with the function blmpSDPD
.
The function SDPMm
provides estimation of a SAR and SDM
model with the Lee-Yu transformation approach (Yu, De Jong, and Lee (2008), Lee and Yu (2010b), Lee
and Yu (2010a)).
The Cigar data set (Baltagi and Levin
(1992)) from the plm
package (Croissant and Millo (2008)) will be used for
describing the use of the blmpSDPD
and SDPMm
functions. It contains panel data of cigarette consumption in 46 states
in the USA over the period from 1963 to 1992. The binary contiguity
matrix of the 46 states is included in the SDPMmod
package.
data("Cigar",package = "plm")
head(Cigar)
#> state year price pop pop16 cpi ndi sales pimin
#> 1 1 63 28.6 3383 2236.5 30.6 1558.305 93.9 26.1
#> 2 1 64 29.8 3431 2276.7 31.0 1684.073 95.4 27.5
#> 3 1 65 29.8 3486 2327.5 31.5 1809.842 98.5 28.9
#> 4 1 66 31.5 3524 2369.7 32.4 1915.160 96.4 29.5
#> 5 1 67 31.6 3533 2393.7 33.4 2023.546 95.5 29.6
#> 6 1 68 35.6 3522 2405.2 34.8 2202.486 88.4 32.0
<- Cigar
data1$logc<-log(data1$sales)
data1$logp<-log(data1$price/data1$cpi)
data1$logy<-log(data1$ndi/data1$cpi)
data1$lpm<-log(data1$pimin/data1$cpi)
data1
data("usa46",package="SDPDmod") ## binary contiguity matrix of 46 USA states
str(usa46)
#> num [1:46, 1:46] 0 0 0 0 0 0 0 1 1 0 ...
#> - attr(*, "dimnames")=List of 2
#> ..$ : chr [1:46] "Alabama" "Arizona" "Arkansas" "California" ...
#> ..$ : chr [1:46] "Alabama" "Arizona" "Arkansas" "California" ...
<- rownor(usa46) ## row-normalization
W isrownor(W) ## check if W is row-normalized
#> [1] TRUE
If the option dynamic
is not set, then the default model
is static. For spatial fixed effects effect
should be set
to “individual”.
<-blmpSDPD(formula = logc ~ logp+logy, data = data1, W = W,
res1index = c("state","year"),
model = list("ols","sar","sdm","sem","sdem","slx"),
effect = "individual")
res1#> $lmarginal
#> ols sar sdm sem sdem slx
#> 1 884.7551 940.015 1049.968 995.7236 1044.108 930.6248
#>
#> $probs
#> ols sar sdm sem sdem slx
#> 1 1.769006e-72 1.765197e-48 0.9971557 2.758904e-24 0.002844276 1.474619e-52
<-blmpSDPD(formula = logc ~ logp+logy, data = data1, W = W,
res2index = c("state","year"),
model = list("ols","sar","sdm","sem","sdem","slx"),
effect = "time")
The default prior is uniform. With prior="beta"
the beta
prior is used.
<-blmpSDPD(formula = logc ~ logp+logy, data = data1, W = W,
res3index = c("state","year"),
model = list("sar","sdm","sem","sdem"),
effect = "twoways",
prior = "beta")
<-blmpSDPD(formula = logc ~ logp+logy, data = data1, W = W,
res4index = c("state","year"),
model = list("sar","sdm","sem","sdem","slx"),
effect = "twoways",
ldet = "mc", ## log-determinant calculated with mcmc procedure
dynamic = TRUE,
prior = "uniform")
<- plm::pdata.frame(data1, index=c('state', 'year'))
d2 $llogc<-plm::lag(d2$logc) ## add lagged variable
d2<-d2[which(!is.na(d2$llogc)),]
data2rownames(data2)<-1:nrow(data2)
<-which(colnames(data2)=="llogc"); kk
kk#> [1] 14
<-blmpSDPD(formula = logc ~ logp+logy, data = data2, W = W,
res5index = c("state","year"),
model = list("sar","sdm","sem","sdem"),
effect = "individual",
ldet = "full",
dynamic = TRUE,
tlaginfo = list(ind=kk),
prior = "beta")
<-SDPDm(formula = logc ~ logp+logy, data = data1, W = W,
mod1index = c("state","year"),
model = "sar",
effect = "individual")
summary(mod1)
#> sar panel model with individual fixed effects
#>
#> Call:
#> SDPDm(formula = logc ~ logp + logy, data = data1, W = W, index = c("state",
#> "year"), model = "sar", effect = "individual")
#>
#> Spatial autoregressive coefficient:
#> Estimate Std. Error t-value Pr(>|t|)
#> rho 0.309172 0.028489 10.852 < 2.2e-16 ***
#>
#> Coefficients:
#> Estimate Std. Error t-value Pr(>|t|)
#> logp -0.52427019 0.02552421 -20.5401 <2e-16 ***
#> logy 0.00024188 0.01518814 0.0159 0.9873
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$rsqr
mod1#> [1] 0.8681613
dynamic
should be set to TRUE as well as the option
stl
in tlaginfo
for inclusion of time lag. For
space-time lag effect the option stl
in
tlaginfo
should also be set to TRUE
<-SDPDm(formula = logc ~ logp+logy, data = data1, W = W,
mod2index = c("state","year"),
model = "sar",
effect = "individual",
dynamic = T,
tlaginfo = list(ind = NULL, tl = T, stl = T))
summary(mod2)
#> sar dynamic panel model with individual fixed effects
#>
#> Call:
#> SDPDm(formula = logc ~ logp + logy, data = data1, W = W, index = c("state",
#> "year"), model = "sar", effect = "individual", dynamic = T,
#> tlaginfo = list(ind = NULL, tl = T, stl = T))
#>
#> Spatial autoregressive coefficient:
#> Estimate Std. Error t-value Pr(>|t|)
#> rho 0.305224 0.031721 9.622 < 2.2e-16 ***
#>
#> Coefficients:
#> Estimate Std. Error t-value Pr(>|t|)
#> logc(t-1) 0.8694130 0.0130195 66.7777 < 2.2e-16 ***
#> W*logc(t-1) -0.2780233 0.0339997 -8.1772 2.905e-16 ***
#> logp -0.1139409 0.0139383 -8.1746 2.967e-16 ***
#> logy -0.0204722 0.0079888 -2.5626 0.01039 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
<-SDPDm(formula = logc ~ logp+logy, data = data1, W = W,
mod3index = c("state","year"),
model = "sar",
effect = "individual",
LYtrans = T,
dynamic = T,
tlaginfo = list(ind = NULL, tl = T, stl = T))
summary(mod3)
#> sar dynamic panel model with individual fixed effects
#>
#> Call:
#> SDPDm(formula = logc ~ logp + logy, data = data1, W = W, index = c("state",
#> "year"), model = "sar", effect = "individual", dynamic = T,
#> tlaginfo = list(ind = NULL, tl = T, stl = T), LYtrans = T)
#>
#> Spatial autoregressive coefficient:
#> Estimate Std. Error t-value Pr(>|t|)
#> rho 0.310607 0.031832 9.7578 < 2.2e-16 ***
#>
#> Coefficients:
#> Estimate Std. Error t-value Pr(>|t|)
#> logc(t-1) 0.9284102 0.0132288 70.1811 < 2.2e-16 ***
#> W*logc(t-1) -0.3017744 0.0354417 -8.5147 2.905e-16 ***
#> logp -0.0858923 0.0139171 -6.1717 2.967e-16 ***
#> logy -0.0215763 0.0081247 -2.6556 0.01039 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
<-SDPDm(formula = logc ~ logp+logy, data = data2, W = W,
mod4index = c("state","year"),
model = "sar",
effect = "individual",
LYtrans = T,
dynamic = T,
tlaginfo = list(ind = kk, tl = T, stl = F))
summary(mod4)
#> sar dynamic panel model with individual fixed effects
#>
#> Call:
#> SDPDm(formula = logc ~ logp + logy, data = data2, W = W, index = c("state",
#> "year"), model = "sar", effect = "individual", dynamic = T,
#> tlaginfo = list(ind = kk, tl = T, stl = F), LYtrans = T)
#>
#> Spatial autoregressive coefficient:
#> Estimate Std. Error t-value Pr(>|t|)
#> rho 0.097235 0.016611 5.8535 1.625e-08 ***
#>
#> Coefficients:
#> Estimate Std. Error t-value Pr(>|t|)
#> logc(t-1) 0.9204780 0.0136652 67.3590 < 2.2e-16 ***
#> logp -0.0498549 0.0141149 -3.5321 1.547e-10 ***
#> logy -0.0316792 0.0083823 -3.7793 0.0002411 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
<-SDPDm(formula = logc ~ logp+logy, data = data1, W = W,
mod5index = c("state","year"),
model = "sdm",
effect = "twoways",
LYtrans = T,
dynamic = T,
tlaginfo = list(ind = NULL, tl = T, stl = T))
summary(mod5)
#> sdm dynamic panel model with twoways fixed effects
#>
#> Call:
#> SDPDm(formula = logc ~ logp + logy, data = data1, W = W, index = c("state",
#> "year"), model = "sdm", effect = "twoways", dynamic = T,
#> tlaginfo = list(ind = NULL, tl = T, stl = T), LYtrans = T)
#>
#> Spatial autoregressive coefficient:
#> Estimate Std. Error t-value Pr(>|t|)
#> rho 0.167440 0.032308 5.1827 0.9097
#>
#> Coefficients:
#> Estimate Std. Error t-value Pr(>|t|)
#> logc(t-1) 0.954193 0.013121 72.722 < 2.2e-16 ***
#> W*logc(t-1) 0.087075 0.040614 2.144 0.9878
#> logp -0.294357 0.025525 -11.532 2.073e-12 ***
#> logy -0.835242 0.032583 -25.634 0.5045
#> W*logp 1.624827 0.044656 36.386 0.6549
#> W*logy 1.608739 0.040520 39.702 0.8993
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
<- impactsSDPDm(mod5)
imp
imp#> $DIRECTst.tab
#> Estimate Std. Error t-value Pr(>|t|)
#> logp -0.2228593 0.02359099 -9.446797 3.493582e-21
#> logy -0.7732413 0.02967081 -26.060672 1.018337e-149
#>
#> $INDIRECTst.tab
#> Estimate Std. Error t-value Pr(>|t|)
#> logp 1.832629 0.05857001 31.28954 6.474728e-215
#> logy 1.706220 0.04958072 34.41298 1.612816e-259
#>
#> $TOTALst.tab
#> Estimate Std. Error t-value Pr(>|t|)
#> logp 1.6097696 0.0654223 24.60583 1.094266e-133
#> logy 0.9329789 0.0445271 20.95306 1.759621e-97
#>
#> $DIRECTlt.tab
#> Estimate Std. Error t-value Pr(>|t|)
#> logp -5.436949 4.26945 -1.27345429 0.2028569
#> logy 2.629254 88.67432 0.02965068 0.9763456
#>
#> $INDIRECTlt.tab
#> Estimate Std. Error t-value Pr(>|t|)
#> logp -0.9756209 4.276854 -0.22811649 0.8195557
#> logy -6.3560624 88.724988 -0.07163779 0.9428902
#>
#> $TOTALlt.tab
#> Estimate Std. Error t-value Pr(>|t|)
#> logp -6.412570 0.5805321 -11.046021 2.291485e-28
#> logy -3.726808 0.4491265 -8.297903 1.059651e-16
#>
#> attr(,"class")
#> [1] "impactsSDPDm"