Sampling from Multi variate Normal with time varying Covariance Matrix in a way similar to matlab

Ooh, yeah, I’m wrong then. You’ll have to keep the loop unless there are math games to be played.

Ignoring the ts again, if you wanted to try sampling the latent variables here,

parameters {
  vector[N] Bz;
}

Bz ~ multi_normal(rep_vector(0, N), B' * diag(A) * B)
y ~ normal(Bz, seps2)

could be done with a non-centered parameterization as (assuming B is NxN matrix) as:

parameters {
  vector[N] z;
}

model {
  vector[N] Bz;
  z ~ normal(0, 1);
  Bz = B' * diag(sqrt(A)) * z;
  y ~ normal(Bz, seps2)
}

I think (just going from: https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution)? But maybe sampling those latents is as painful as doing the decompositions all the time. I dunno.