Hidden Markov model with missing values

hmm_missing.stan

data {
  int<lower = 0> N_ind;
  int<lower = 0> N_t;
  int<lower = 0, upper = 2> Y[N_ind, N_t]; // 0: missing value
}

parameters {
  real<lower = 0, upper = 1> phi;
  real<lower = 0, upper = 1> psi;
  real<lower = 0, upper = 1> p;
}

transformed parameters {
  simplex[2] ps[2];
  simplex[2] po[2];

  ps[1, 1] = phi;
  ps[1, 2] = 1 - phi;
  ps[2, 1] = 1 - psi;
  ps[2, 2] = psi;

  po[1, 1] = p;
  po[1, 2] = 1 - p;
  po[2, 1] = 0;
  po[2, 2] = 1;
}

model {
  real acc[2];
  vector[2] gamma[N_t];

  for (i in 1:N_ind) {
    for (k in 1:2)
      gamma[1, k] = (k == Y[i, 1]);

    for (t in 2:N_t) {
      for (k in 1:2) {
        for (j in 1:2) {
          if (Y[i, t] == 0) {
            acc[j] = gamma[t - 1, j] * ps[j, k];
          } else {
            acc[j] = gamma[t - 1, j] * ps[j, k] * po[k, Y[i, t]];
          }
        }
        gamma[t, k] = sum(acc);
      }
    }
    target += log(sum(gamma[N_t]));
  }
}