Generating Dynamic Linear Regression Coefficients

Hello, I am quite new to Stan. I have drafted a Kalman filter-based dynamic linear regression model as below. After running the model, I can obtain the posterior distributions for those model parameters, but I am struggling on how to obtain the regression coefficients based on those posterior distributions of model parameters. I am considering about adding a generated quantities block, but failed to do so. Any help would be appreciatived. Thank you so much.

data {
  // number of observations; In our case, n = 18 years
  int n;
  // number of variable == 413; In our case, p = 413 as we have 413 block groups
  int p;
  // number of states; In our case, we have 5 predictors for the housing price, thus m = 5
  int m;
  // size of system error; In our case, the size of system error should equal to 5
  int r;
  // observed data; In our case, the observed data y is our time-series housing price data
  vector[p] y[n];
  // observation eq; In our case, the observation eq is our predictors
  matrix[p, m] ZZ[n];
  // system eq; T is the evolution matrix. 
  //matrix[m, m] T; In our case, T is unknown, and will be estimated as parameter
  matrix[m, r] R;
  // initial values; In our case, the a_1 and P_a is unknow parameter
  //vector[m] a_1;
  //cov_matrix[m] P_1;
  // measurement error
  // vector<lower=0.0>[p] h;
  // vector<lower=0.0>[r] q; 
  int<lower = 0, upper = 1> run_estimation;
}
transformed data { 
  real log_2p = log(2 * pi());
}
parameters {
  vector<lower=0.0>[p] h;
  vector<lower=0.0>[r] q;
  vector[m] t;
  vector[m] a_1;
  corr_matrix[m] Omega;
  vector<lower=0.0>[m] tau;
}
transformed parameters {
 // Unless you explicitly want a matrix in the return / don't wanna get a vector back
 // matrix[m, m] T = diag_matrix(t);
  // Unless you want this in the return
  //cov_matrix[m] P_1 = quad_form_diag(Omega, tau);  
}
model {
  vector[m] a[n + 1];
  matrix[m, m] P[n + 1];
  real llik_obs[n];
  real llik;

  //Set priors
  tau ~ normal(0, 10);
  Omega ~ lkj_corr(2);

  h ~ normal(0, 10);
  q ~ normal(0, 10);
  t ~ normal(0, 10);
  a_1 ~ normal(0, 10);
  // 1st observation
  a[1] = a_1;
  P[1] = quad_form_diag(Omega, tau);

  for (i in 1:n) {
    vector[p] v = y[i] - ZZ[i] * a[i];
    matrix[p, p] F = ZZ[i] * P[i] * ZZ[i]'+ diag_matrix(h);
    matrix[p, p] Finv = inverse(F);
    matrix[m, p] K = diag_pre_multiply(t, P[i]) * ZZ[i]' * Finv;
    // // manual update of multivariate normal
    llik_obs[i] = -0.5 * (p * log_2p + log(determinant(F)) + v' * Finv * v);
  
    a[i + 1] = t .* a[i] + K * v;
    //a[i + 1] = T * a[i] + K * v;
    P[i + 1] = diag_pre_multiply(t, P[i]) * (diag_matrix(t) - K * ZZ[i])' + diag_post_multiply(R, q) * R';
  }
  llik = sum(llik_obs);
  
  if(run_estimation == 1){
    target += llik; 
  }
}


1 Like

Afternoon,

Could you post your attempt at the generate quantity block? Thanks!

Thank you. The generated quantities block is as below:

generated quantities {
    vector[m] coef[n];
    coef[1] = a_1;
    for (i in 2:n){
     coef[i] = multi_normal_rng(diag_matrix(t)*coef[i-1], diag_matrix(q));
    }
}
2 Likes