I stumpled upon some really good answers in this forum so I decided to give it a try and hope that I can also enrich the ongoing discussions here.

I try to implement a multivariate stochastic volatility model in Winbugs. Hereby I use the matrix Sigma.epsilon with stochastic parts in the off-diagonal elements (correlations) and deterministic values in the diagonals.

However, when investigating the draws of Sigma.epsilon after sampling, Winbugs returns only matrices with 0's in the diagonal.

I do not see where the mistake could be, every comment is appreciated

Code:

```
model{
for (t in 1:T)
{
Omega[t,1,1] <- exp(h[1,t]);
Omega[t,2,2] <- exp(h[2,t]);
Omega[t,1,2] <- 0;
Omega[t,2,1] <- 0;
Sigma.epsilon[t,1,1] <- 1;
Sigma.epsilon[t,2,2] <- 1;
Sigma.epsilon[t,1,2] <- rho.e;
Sigma.epsilon[t,2,1] <- rho.e;
for (n in 1:N){
for (j in 1:N){
Var.y[t,n,j] <- Omega[t,n,j]*Sigma.epsilon[t,n,j]*Omega[t,j,n];
}}
Prec.y[t,1:N,1:N]<-inverse(Var.y[t,,])
y[t,1:2] ~ dmnorm(mu.y[1:2],Prec.y[t,,]);
}
for (n in 1:N)
{
h.st[1,n] <- mu.h[n]
for (t in 2:T)
{
h.st[t,n] ~dnorm(mu.h[n],1);
h.mu[t,n] <- mu.h[n] + ph[n]*(h.st[t-1,n]-mu.h[n]);
}
}
for (t in 1:T)
{
h[1,t] <- sig.u[1]*h.st[t,1];
h[2,t] <- sig.u[2]*rho.u*h.st[t,1]+sig.u[2]*sqrt(1-rho.u*rho.u)*h.st[t,2];
}
# priors
for (n in 1:N)
{
invsig2u[n] ~dgamma(2.5,0.025);
sig.u[n] <- sqrt(1/invsig2u[n]);
phstar[n] ~dbeta(20,1.5);
ph[n] <- 2*phstar[n] -1;
mu.h[n] ~dnorm(0,0.04);
}
rho.e ~dunif(-1,1);
rho.u ~dunif(0,1);
mu.y[1] <- 0;
mu.y[2] <- 0;
}
```