# Fixed effects spatial panel data models with time-varying spatial autocorrelation parameters

Hi there,

I want to estimate a spatial lag panel data model with time-varying spatial autocorrelation parameters.

y_{it}=\rho_tWy_{it}+x_{it}\beta +\epsilon , \epsilon \sim N(0, \sigma^2)

I use the following rstan code, but I did not get the correct estimates. The fact is that \rho_t is always under-estimated.
I wonder Is there any way to improve my code? Or is there anything wrong in my code.

Many thanks
C.LING

functions {
matrix kronecker_prod(matrix A, matrix B) {
matrix[rows(A) * rows(B), cols(A) * cols(B)] C;
int m;
int n;
int p;
int q;
m = rows(A);
n = cols(A);
p = rows(B);
q = cols(B);
for (i in 1:m)
for (j in 1:n)
for (k in 1:p)
for (l in 1:q)
C[p*(i-1)+k,q*(j-1)+l] <- A[i,j]*B[k,l];
return C;
}
/* normal log-pdf for spatially lagged responses
* Args:
*   y: the response vector
*   mu: mean parameter vector
*   sigma: residual standard deviation
*   rho: a vector of autoregressive parameters
*   W: spatial weight matrix
* Returns:
*   a scalar to be added to the log posterior
*/
real normal_lagsar_lpdf(vector y, vector mu, real sigma,
vector rho, data matrix W) {
int NT = rows(y);
int N = rows(W);
int T = NT%/%N;
real inv_sigma2 = inv_square(sigma);
matrix[NT, NT] IW_tilde = add_diag(-kronecker_prod(diag_matrix(rho), W), 1);
vector[NT] half_pred;
real log_det;
half_pred = IW_tilde * y - mu;
log_det = log_determinant(IW_tilde);
return  0.5 * NT * log(inv_sigma2) + log_det -
0.5 * dot_self(half_pred) * inv_sigma2;
}
}

data {
int<lower=1> NT;  // total number of observations
int<lower=1> T;  // total number of periods
int<lower=1> N;  // total number of observations in a cross-section
vector[NT] Y;  // response variable
int<lower=1> K;  // number of population-level effects
matrix[NT, K] X;  // population-level design matrix
matrix[N, N] Msar;  // spatial weight matrix
int prior_only;  // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[NT, Kc] Xc;  // centered version of X without an intercept
vector[Kc] means_X;  // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b;  // population-level effects
real Intercept;  // temporary intercept for centered predictors
vector<lower=-1,upper=1>[T] lagsar;  // lag-SAR correlation parameters
real<lower=0> sigma;  // dispersion parameter
}
transformed parameters {
real lprior = 0;  // prior contributions to the log posterior
lprior += student_t_lpdf(Intercept | 3, 34, 21.6);
lprior += student_t_lpdf(sigma | 3, 0, 21.6)
- 1 * student_t_lccdf(0 | 3, 0, 21.6);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[NT] mu = rep_vector(0.0, NT);
mu += Intercept + Xc * b;
target += normal_lagsar_lpdf(Y | mu, sigma, lagsar, Msar);
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}