Hello, I wanna use state space model(ssm) to estimate the parameters. but my sampling is too slow. My rstan code:
functions{
int symmat_size(int n) {
int sz;
// This calculates it iteratively because Stan gives a warning
// with integer division.
sz = 0;
for (i in 1:n) {
sz = sz + i;
}
return sz;
}
matrix to_symmetric_matrix(matrix x) {
return 0.5 * (x + x ');
}
vector symmat_to_vector(matrix x) {
vector[symmat_size(rows(x))] v;
int k;
k = 1;
// if x is m x n symmetric, then this will return
// only parts of an m x m matrix.
for (j in 1:rows(x)) {
for (i in 1:j) {
v[k] = x[i, j];
k = k + 1;
}
}
return v;
}
vector ssm_filter_update_a(vector a, vector c, matrix T, vector v, matrix K) {
vector[num_elements(a)] a_new;
a_new = T * a + K * v + c;
return a_new;
}
matrix ssm_filter_update_P(matrix P, matrix Z, matrix T,
matrix RQR, matrix K) {
matrix[rows(P), cols(P)] P_new;
P_new = to_symmetric_matrix(T * P * (T - K * Z)' + RQR);
return P_new;
}
vector ssm_filter_update_v(vector y, vector a, vector d, matrix Z) {
vector[num_elements(y)] v;
v = y - Z * a - d;
return v;
}
matrix ssm_filter_update_F(matrix P, matrix Z, matrix H) {
matrix[rows(H), cols(H)] F;
F = quad_form(P, Z') + H;
return F;
}
matrix ssm_filter_update_Finv(matrix P, matrix Z, matrix H) {
matrix[rows(H), cols(H)] Finv;
Finv = inverse(ssm_filter_update_F(P, Z, H));
return Finv;
}
matrix ssm_filter_update_K(matrix P, matrix Z, matrix T, matrix Finv) {
matrix[cols(Z), rows(Z)] K;
K = T * P * Z' * Finv;
return K;
}
real ssm_filter_update_ll(vector v, matrix Finv) {
real ll;
int p;
p = num_elements(v);
// det(A^{-1}) = 1 / det(A) -> log det(A^{-1}) = - log det(A)
ll = (- 0.5 *
(p * log(2 * pi())
- log_determinant(Finv)
+ quad_form(Finv, v)
));
return ll;
}
real ssm_filter_get_loglik(vector x, int m, int p) {
real y;
y = x[1];
return y;
}
real ssm_constant_lpdf(vector[] y,
vector d, matrix Z, matrix H,
vector c, matrix T, matrix Q,
vector a1, matrix P1) {
real ll;
int n;
int m;
int p;
n = size(y); // number of obs
m = cols(Z);
p = rows(Z);
{
vector[n] ll_obs;
vector[m] a;
matrix[m, m] P;
vector[p] v;
matrix[p, p] Finv;
matrix[m, p] K;
a = a1;
P = P1;
for (t in 1:n) {
v = ssm_filter_update_v(y[t], a, d, Z);
Finv = ssm_filter_update_Finv(P, Z, H);
K = ssm_filter_update_K(P, Z, T, Finv);
ll_obs[t] = ssm_filter_update_ll(v, Finv);
if (t < n) {
a = ssm_filter_update_a(a, c, T, v, K);
// check for convergence
// should only check for convergence if there are no missing values
P = ssm_filter_update_P(P, Z, T, Q, K);
}
}
ll = sum(ll_obs);
}
return ll;
}
}
data {
int<lower=1> WT; //the number of sample
int<lower=1> WK; //number of K years
//real y1[T]; //data1
vector[11] y[WT]; //data2
//matrix[T,K] y3; // data3
real h;
vector[WK] tau; // represent K different year
}
parameters {
//parameters for mu
vector<lower=0>[5] kappa;
vector[4] kappaq;
vector[5] mu;
vector<lower=0>[5] sigma;
vector<lower=0>[11] eps;
real cr;
real ca;
real cg;
real gr;
}
transformed parameters{
//matrix[5,5] sigr=diag_matrix(sigma.*sigma);
matrix[5,5] phi;
cov_matrix[5] Q;
// matrix[5,WT] muF;
matrix[5,5] sigF1;
matrix[11,5] Z=rep_matrix(0,11,5);
matrix[5,4] C1a;
vector[WK] Ca;
vector[WK] Cg;
vector[11] C0;
cov_matrix [11] H;
vector[5] muu;
H = diag_matrix(eps);
phi = diag_matrix(exp(-(kappa)*h));
Q = diag_matrix(exp(2*log(sigma)).*(1-exp(-2*kappa*h))./(2*kappa));
//muF[,1] =mu;
sigF1 = diag_matrix(exp(2*log(sigma))./(2*kappa));
for(i in 1:5){
for (j in 1:4){
C1a[i,j] = (1-exp(-kappaq[j]*tau[i]))/(kappaq[j]*tau[i]);
}
}
for (j in 1:WK){
Ca[j]=(cr+ca) - sigma[1]^2 *(1-C1a[j,1])/(2*kappaq[1]^2)+sigma[1]^2 *C1a[j,1]^2*tau[j]/(4*kappaq[1])-
sigma[2]^2 *(1-C1a[j,2])/(2*kappaq[2]^2)+sigma[2]^2 *C1a[j,2]^2*tau[j]/(4*kappaq[2])-
sigma[3]^2 *(1-C1a[j,3])/(2*kappaq[3]^2)+sigma[3]^2 *C1a[j,3]^2*tau[j]/(4*kappaq[3]);
}
for(j in 1:WK){
Cg[j]=(cr*gr+cg) - gr^2*sigma[1]^2 *(1-C1a[j,1])/(2*kappaq[1]^2)+gr^2*sigma[1]^2 *C1a[j,1]^2*tau[j]/(4*kappaq[1])-
gr^2*sigma[2]^2 *(1-C1a[j,2])/(2*kappaq[2]^2)+gr^2*sigma[2]^2 *C1a[j,2]^2*tau[j]/(4*kappaq[2])-
sigma[4]^2 *(1-C1a[j,4])/(2*kappaq[4]^2)+sigma[4]^2 *C1a[j,4]^2*tau[j]/(4*kappaq[4]);
}
C0[1]=cr;
C0[2:6]=Ca;
C0[7:11]=Cg;
Z[1,1] = 1;
Z[1,2] = 1;
Z[1,5] = 1;
for (i in 2:6){
for(j in 1:3){
Z[i,j] = C1a[i-1,j];
}
}
for (i in 1:5){
Z[i+6,1] = C1a[i,1]*gr;
Z[i+6,2] = C1a[i,2]*gr;
Z[i+6,4] = C1a[i,4];
}
muu=mu-phi*mu;
}
model {
vector[ssm_filter_size(cols(Z), rows(Z))] res1[size(y)];
kappa~normal(0,1);
kappaq~normal(0,5);
mu ~normal(0,0.02);
sigma ~normal(0,1);
eps ~ normal(0,0.02);
cr ~ normal(0,0.02);
ca ~ normal(0,0.02);
cg ~ normal(0,0.02);
gr ~ normal(0,1);
target += ssm_constant_lpdf(y|C0,Z,H,muu, phi, Q, mu, sigF1);
}
Some suggestions?
very thanks!