Hi all,

I am using rstan to fit a bivariate multilevel VAR model. Since my model has two random intercepts, two random autoregressive parameters, and two random cross-regression parameters, my random effect covariance matrix is a 6 dimensional matrix, following a multivariate normal distribution.

I understand that the preferred way of parameterizing the covariance matrix in Stan is to use Cholesky decomposition since it is more efficient and numerical stable. I did obtain good recovery based on simulated data, following this parameterization. However, when I defined the covariance matrix as â€ścov_matrix[6] bcovâ€ť and passed it to â€śmulti_normal(mu, bcov)â€ť, and assigned an inverse Wishart prior to bcov, the random effect variances were **overestimated**. I checked the mean, median, and mode of the posterior samples, and they were similar, so the posterior samples were not skewed. The model has converged and the effective sample sizes were large (e.g., about 10k).

So I am just curious about why it is the case. I can provide more details and/or code if necessary.

Thanks in advance.