I tried to apply the state space model using Filtering and Smoothing as the original codes are written by Jeff Arnold from his Githubs.
my file is named as “ssmfunc.stan” where the code contains three functions namely 1) ssm_lpdf, 2) ssm_filter_size, 3) ssm_simsmo_states_rng, 4) ssm_filter
real ssm_lpdf(vector[] y,
vector[] d, matrix[] Z, matrix[] H,
vector[] c, matrix[] T, matrix[] R, matrix[] Q,
vector a1, matrix P1) {
real ll;
int n;
int m;
int p;
int q;
n = size(y);
m = dims(Z)[3];
p = dims(Z)[2];
q = dims(Q)[2];
{
vector[p] d_t;
matrix[p, m] Z_t;
matrix[p, p] H_t;
vector[m] c_t;
matrix[m, m] T_t;
matrix[m, q] R_t;
matrix[q, q] Q_t;
matrix[m, m] RQR;
vector[n] ll_obs;
vector[m] a;
matrix[m, m] P;
vector[p] v;
matrix[p, p] Finv;
matrix[m, p] K;
d_t = d[1];
Z_t = Z[1];
H_t = H[1];
c_t = c[1];
T_t = T[1];
R_t = R[1];
Q_t = Q[1];
RQR = quad_form_sym(Q_t, R_t ');
a = a1;
P = P1;
for (t in 1:n) {
if (t > 1) {
if (size(d) > 1) {
d_t = d[t];
}
if (size(Z) > 1) {
Z_t = Z[t];
}
if (size(H) > 1) {
H_t = H[t];
}
if (t < n) {
if (size(c) > 1) {
c_t = c[t];
}
if (size(T) > 1) {
T_t = T[t];
}
if (size(R) > 1) {
R_t = R[t];
}
if (size(Q) > 1) {
Q_t = Q[t];
}
if (size(R) > 1 || size(Q) > 1) {
RQR = quad_form_sym(Q_t, R_t ');
}
}
}
v = ssm_update_v(y[t], a, d_t, Z_t);
Finv = ssm_update_Finv(P, Z_t, H_t);
K = ssm_update_K(P, T_t, Z_t, Finv);
ll_obs[t] = ssm_update_loglik(v, Finv);
if (t < n) {
a = ssm_update_a(a, c_t, T_t, v, K);
P = ssm_update_P(P, Z_t, T_t, RQR, K);
}
}
ll = sum(ll_obs);
}
return ll;
}
vector[] ssm_filter(vector[] y,
vector[] d, matrix[] Z, matrix[] H,
vector[] c, matrix[] T, matrix[] R, matrix[] Q,
vector a1, matrix P1) {
vector[ssm_filter_size(dims(Z)[3], dims(Z)[2])] res[size(y)];
int q;
int n;
int p;
int m;
n = size(y);
m = dims(Z)[3];
p = dims(Z)[2];
q = dims(Q)[2];
{
vector[p] d_t;
matrix[p, m] Z_t;
matrix[p, p] H_t;
vector[m] c_t;
matrix[m, m] T_t;
matrix[m, q] R_t;
matrix[q, q] Q_t;
matrix[m, m] RQR;
vector[m] a;
matrix[m, m] P;
vector[p] v;
matrix[p, p] Finv;
matrix[m, p] K;
real ll;
int idx[6, 3];
idx = ssm_filter_idx(m, p);
d_t = d[1];
Z_t = Z[1];
H_t = H[1];
c_t = c[1];
T_t = T[1];
R_t = R[1];
Q_t = Q[1];
RQR = quad_form_sym(Q_t, R_t ');
a = a1;
P = P1;
for (t in 1:n) {
if (t > 1) {
if (size(d) > 1) {
d_t = d[t];
}
if (size(Z) > 1) {
Z_t = Z[t];
}
if (size(H) > 1) {
H_t = H[t];
}
if (size(c) > 1) {
c_t = c[t];
}
if (size(T) > 1) {
T_t = T[t];
}
if (size(R) > 1) {
R_t = R[t];
}
if (size(Q) > 1) {
Q_t = Q[t];
}
if (size(R) > 1 || size(Q) > 1) {
RQR = quad_form_sym(Q_t, R_t ');
}
}
v = ssm_update_v(y[t], a, d_t, Z_t);
Finv = ssm_update_Finv(P, Z_t, H_t);
K = ssm_update_K(P, Z_t, T_t, Finv);
ll = ssm_update_loglik(v, Finv);
res[t, 1] = ll;
res[t, idx[2, 2]:idx[2, 3]] = v;
res[t, idx[3, 2]:idx[3, 3]] = symmat_to_vector(Finv);
res[t, idx[4, 2]:idx[4, 3]] = to_vector(K);
res[t, idx[5, 2]:idx[5, 3]] = a;
res[t, idx[6, 2]:idx[6, 3]] = symmat_to_vector(P);
if (t < n) {
a = ssm_update_a(a, c_t, T_t, v, K);
P = ssm_update_P(P, Z_t, T_t, RQR, K);
}
}
}
return res;
}
vector[] ssm_simsmo_states_rng(vector[] filter,
vector[] d, matrix[] Z, matrix[] H,
vector[] c, matrix[] T, matrix[] R, matrix[] Q,
vector a1, matrix P1) {
vector[dims(Z)[3]] draws[size(filter)];
int n;
int p;
int m;
int q;
n = size(filter);
m = dims(Z)[3];
p = dims(Z)[2];
q = dims(Q)[2];
{
vector[ssm_filter_size(m, p)] filter_plus[n];
vector[ssm_sim_size(m, p, q)] sims[n];
vector[p] y[n];
vector[m] alpha_hat_plus[n];
vector[m] alpha_hat[n];
alpha_hat = ssm_smooth_states_mean(filter, Z, c, T, R, Q);
sims = ssm_sim_rng(n, d, Z, H, c, T, R, Q, a1, P1);
for (i in 1:n) {
y[i] = ssm_sim_get_y(sims[i], m, p, q);
}
filter_plus = ssm_filter(y, d, Z, H, c, T, R, Q, a1, P1);
alpha_hat_plus = ssm_smooth_states_mean(filter_plus, Z, c, T, R, Q);
for (i in 1:n) {
draws[i] = (ssm_sim_get_a(sims[i], m, p, q)
- alpha_hat_plus[i]
+ alpha_hat[i]);
}
}
return draws;
}
int ssm_filter_size(int m, int p) {
int sz;
int idx[6, 3];
idx = ssm_filter_idx(m, p);
sz = idx[6, 3];
return sz;
}
and Here is the code to for compiling (notice that i havn’t edit any single letter from the original Jeff Arnold’s code:
functions {
#include /ssmfunc.stan
}
data {
int<lower = 1> n;
vector[1] y[n];
int<lower = 1> k;
vector[k] x[n];
vector<lower = 0.0>[1] a1;
cov_matrix[1] P1;
real<lower = 0.0> y_scale;
}
transformed data {
// system matrices
matrix[1, 1] T[1];
matrix[1, 1] Z[1];
matrix[1, 1] R[1];
vector[1] c[1];
int m;
int p;
int q;
int filter_sz;
m = 1;
p = 1;
q = 1;
T[1, 1, 1] = 1.0;
Z[1, 1, 1] = 1.0;
R[1, 1, 1] = 1.0;
c[1, 1] = 0.0;
filter_sz = ssm_filter_size(m, p);
}
parameters {
real<lower = 0.0> sigma_eta;
real<lower = 0.0> sigma_epsilon;
vector[k] beta;
}
transformed parameters {
matrix[1, 1] H[1];
matrix[1, 1] Q[1];
vector[1] d[n];
H = rep_array(rep_matrix(pow(sigma_epsilon, 2), 1, 1), 1);
Q = rep_array(rep_matrix(pow(sigma_eta * sigma_epsilon, 2), 1, 1), 1);
for (t in 1:n) {
d[t, 1] = dot_product(beta, x[t]);
}
}
model {
target += ssm_lpdf(y | d, Z, H, c, T, R, Q, a1, P1);
sigma_epsilon ~ cauchy(0.0, y_scale);
sigma_eta ~ cauchy(0.0, 1.0);
}
generated quantities {
vector[filter_sz] filtered[n];
vector[1] alpha[n];
vector[1] mu[n];
// filtering
filtered = ssm_filter(y, d, Z, H, c, T, R, Q, a1, P1);
// sampling states
alpha = ssm_simsmo_states_rng(filtered, d, Z, H, c, T, R, Q, a1, P1);
// sample E(y_t | y_1, ..., y_n)
for (t in 1:n) {
mu[t] = d[t] + Z[1] * alpha[t];
}
}
Here is the error i got:
I found the topic Stanc could not find include file - #5 by ron tried it but it did not work.
Any help would be really appreciated.
Best,
PATT