I'm asking about the estimation of direct, indirect, and total effects of spatial durbin panel data modeling. As far as I know, there is no direct package of doing the estimations. So I generate the spatial lag of independent variables myself. Here is the simple version of my code
# Generate the spatial lags:
dataS$JC_2dg_s <- slag(dataS$JC_2dg, listw = cont.listw)
dataS$JD_2dg_s <- slag(dataS$JD_2dg, listw = cont.listw)
# Define the equation and make the estimation
eq0_1 <- TFR_2 ~ JC_2dg + JD_2dg + JC_2dg_s + JD_2dg_s +
FLFP2030 + FLFP2030_s + logGDP_r + logGDP_r_s + expat_ratio + expat_ratio_s +
hp1000 + hp_s
model1<-spml(eq0_1, data=dataS, listw=cont.listw, model="within", spatial.error="none",lag=T, effect="twoway", rel.tol=2e-40)
summary(model1)
## But then, how to get the effects? I try to do it manually:
W.c <- as(as_dgRMatrix_listw(cont.listw), "CsparseMatrix")
beta_jc <- model1$coefficients[2]
theta_jc <- model1$coefficients[4]
beta_jd <- model1$coefficients[3]
theta_jd <- model1$coefficients[5]
# Calculate rho * W
rhoW <- model1$spat.coef * W.c
# Identity matrix
I <- Matrix::Diagonal(nrow(W.c))
# Calculate (I - rho * W)^(-1)
inverse_matrix <- solve(I - rhoW)
# Calculate the direct impact of beta and theta using (I - rho * W)^(-1)(beta + W * theta)
mateff_jc <- inverse_matrix *(beta_jc) + inverse_matrix %*% (W.c *theta_jc)
mateff_jd <- inverse_matrix *(beta_jd) + inverse_matrix %*% (W.c *theta_jd)
# Calculate the average direct effect and its standard error
ade_jc <- mean(diag(mateff_jc))
se_ade_jc <- sd(diag(mateff_jc))
t_directjc <- ade_jc/se_ade_jc
p_directjc <- 2 * (1 - pt(abs(t_directjc), df = nrow(W.c) - 1))
ade_jd <- mean(diag(mateff_jd))
se_ade_jd <- sd(diag(mateff_jd))
t_directjd <- ade_jd/se_ade_jd
p_directjd <- 2 * (1 - pt(abs(t_directjd), df = nrow(W.c) - 1))`
But after checking some literature, I don't know if I made the calculation of standard errors correct, as some variance-covariance might be included in the function.