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

**URL:** <https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343>\
**Category:** Modeling\
**Created:** [February 18, 2018, 3:55pm UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343 "2018-02-18T15:55:43Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![ChrisMu](https://avatars.discourse-cdn.com/v4/letter/c/d2c977/32.png) [@ChrisMu](https://discourse.mc-stan.org/u/ChrisMu)\
**Post date:** [February 18, 2018, 3:55pm UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343/1 "2018-02-18T15:55:43Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [February 19, 2018, 1:44am UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343/2 "2018-02-19T01:44:00Z")

</div>

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);
  ...

```

---

<div class="post-metadata">

**Author:** ![ChrisMu](https://avatars.discourse-cdn.com/v4/letter/c/d2c977/32.png) [@ChrisMu](https://discourse.mc-stan.org/u/ChrisMu)\
**Post date:** [February 19, 2018, 8:16am UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343/3 "2018-02-19T08:16:19Z")

</div>

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[(2_k_N - 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)) ;  
}

}

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [February 19, 2018, 6:30pm UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343/4 "2018-02-19T18:30:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisMu](https://avatars.discourse-cdn.com/v4/letter/c/d2c977/32.png) [@ChrisMu](https://discourse.mc-stan.org/u/ChrisMu)\
**Post date:** [February 22, 2018, 3:40pm UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343/5 "2018-02-22T15:40:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [February 22, 2018, 5:28pm UTC](https://discourse.mc-stan.org/t/sampling-from-multi-variate-normal-with-time-varying-covariance-matrix-in-a-way-similar-to-matlab/3343/6 "2018-02-22T17:28:14Z")

</div>

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

Ignoring the `t`s 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](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.