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