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

Modeling
1

Here is what I am currently coding:

model {
for (j in 1:k) {
h_std[j] ~ normal(0,1) ;
}
for (t in 1:TT) {
y[,t] ~ multi_normal(rep_vector(0,N),quad_form_sym(diag_matrix(to_vector(exp(h[,t]))),B’)+diag_matrix(seps2)) ;
}

}

where h_std is an array of k vectors of length TT

vector[TT] h_std[k] ;

and y is a N times TT matrix containing the data:

matrix[N,TT] y ; // matrix of returns

what I am wondering: is there a way to state the sampling argument for y in a vector notation? E.g. in matlab, when I have a three dimensional array Sigma that is of dimensions N x N x TT, i.e. containing lets say TT = 1000 N x N variance covariance matrices, I can write

Y = mvnrnd(zeros(TT,N), Sigma);

which gives me a matrix of dimension TT x N numbers drawn from the multivariate normal distribution with mean zero and the covariance matrix (or more precise the TT covariance matrices) Sigma.

Is there a way to so some thing similar in stan?

2

It helps if you include your whole Stan model.

If you want to use random number generators, you can do that in the generated quantities block, but it’s not fully vectorized yet (probably two more releases). If it’s in the model block, then yes, you can vectorize.

I couldn’t parse the two arguments to quad_form_sym, but if you stick to diagonal covariance matrices, you should just use normal rather than multi_normal. That is,

y ~ multi_normal(mu, diag_matrix(sigma));

is just a much less efficient version of

y ~ normal(mu, sigma);

You can vectorize your sampling statement if you transpose the ordering for y and make it an array of vectors,

data {
  vector[N] y[TT];
  ...
transformed data {
  vector[N] zeros = rep_vector(0, N);
   ...
model {
  y ~ multi_normal(zeros, Sigma);
  ...
3

Here is the complete code. I work with simulated data. The problem is that the covariance matrix is not diagonal. It is defined as B*exp(h_t)*B’ + seps2, where exp(h_t) is a diagonal matrix as is seps2. Multiplying B and B’ makes the whole thing non-diagonal.

data {
int<lower = 0> TT ; // length of the data
int<lower = 0> N ; // number of time series
int<lower = 0> k ; //number of factors
matrix[N,TT] y ; // matrix of returns
}

parameters {
row_vector<lower = -1 , upper = 1>[k] phi; //autoregressive parameter
// row_vector[k] mu ;//constant of the autoregressive process
vector<lower = 0>[k] seta2 ; // variance of the error in the autoregressive process
vector<lower = 0>[N] seps2; // Defining the diagonal elements of Omega
vector[TT] h_std[k] ;
vector[(2kN - k*k - k)/2] b;

}

transformed parameters {
matrix[N,k] B ; // factor loading
vector[TT] h[k] ;
row_vector[k] phistar ;
row_vector[k] mu;
mu = rep_row_vector(0,k);

for (i in 1:k) {
for (j in i:k) {
B[1:i-1,i] = rep_vector(0,i-1);
}
B[i,i] = 1;
}
{
int tmp ;
tmp = 1;
for (i in 1:k) {
for (j in i:N-1) {
B[j+1,i] = b[tmp] ;
tmp = tmp + 1;
}
}
}
phistar = (phi + 1)/2 ;
for (j in 1:k) {
h[j] = sqrt(seta2[j])*h_std[j] ;
h[j,1] = h[j,1]/sqrt(1-phi[j]*phi[j]);
h[j] = h[j] + mu[j] ;
for (t in 2:TT) {
h[j,t] = h[j,t] + phi[j] * (h[j,t-1] - mu[j]);
}
}

}

model {
for (j in 1:k) {
h_std[j] ~ normal(0,1) ;
}
for (t in 1:TT) {
y[,t] ~ multi_normal(rep_vector(0,N),quad_form_sym(diag_matrix(to_vector(exp(h[,t]))),B’)+diag_matrix(seps2)) ;
}

}

4

What @Bob_Carpenter said is right, the code block at the end of his post will vectorize fine with a non-diagonal covariance matrix.

5

Sorry for the late reply. I was trying to implement your suggestions. The problem is, that my Sigma is time varying. So what I have are T values of log volatilities. These T values have to be multiplied with the quadratic form of the matrix B, therefore I do the quad_form_sym(diag_matrix(to_vector(exp(h[,t]))),B’)
plus a constant variance covariance matrix from the idiosynchratic error
diag_matrix(seps2)
within a loop.
defining y as an array of vectors does not change the fact that I have to use a loop over the time periods for the sampling, at least I did not manage to realise that.

6

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.