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

Modeling
, ,
1

Hi there,

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

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);
}