I am trying to predict median housing price at the block group level using a bunch of predictors (e.g., median household income, distance to the CBD, etc). I have 413 block groups, and each block group has an annual median housing price from 1993 to 2010. I try to fit a dynamic linear regression between the housing price and the predictors using kalman filter. The code is shown as below. When running for just several block groups, the speed is okay. But when using all the 413 block groups, it is extremely slow. I am trying to reparamterize this model, but so far havenâ€™t found any effective solution.

```
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
real ZZ[m, p, 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;
}
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 {
matrix[p, p] H;
matrix[r, r] Q;
matrix[m, m] T;
cov_matrix[m] P_1;
H = diag_matrix(h);
Q = diag_matrix(q);
T = diag_matrix(t);
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);
for (i in 1:p)
h[i] ~ normal(0, 10);
for (i in 1:r)
q[i] ~ normal(0, 10);
for (i in 1:m)
t[i] ~ normal(0, 10);
for (i in 1:m)
a_1[i] ~ normal(0, 10);
// 1st observation
a[1] = a_1;
P[1] = P_1;
for (i in 1:n) {
vector[p] v;
matrix[p, p] F;
matrix[p, p] Finv;
matrix[m, p] K;
matrix[m, m] L;
matrix[p, m] Z;
Z = to_matrix(ZZ[, , n])';
v = y[i] - Z * a[i];
F = Z * P[i] * Z' + H;
Finv = inverse(F);
K = T * P[i] * Z' * Finv;
L = T - K * Z;
// // manual update of multivariate normal
llik_obs[i] = -0.5 * (p * log(2 * pi()) + log(determinant(F)) + v' * Finv * v);
a[i + 1] = T * a[i] + K * v;
P[i + 1] = T * P[i] * L' + R * Q * R';
}
llik = sum(llik_obs);
if(run_estimation == 1){
target += llik;
}
}
```