I’m making a multivariate stochastic volatility model with correlation on the returns but not the volatility processes. I’ve already reparameterized the volatility process parameter ‘h’, which dramatically improved the speed and mixing. I’m wondering if it’s possible to also reparameterize the sampling distribution for the ‘returns’ data variable. It’s my understanding that it should be possible to specify returns as being sampled from a standard normal rescaled by the volatility taken from ‘h’ and the Cholesky decomposition of the correlation matrix, but I’m not sure how to implement that. Any other tips for efficiency are appreciated as well. Thanks!

```
data {
int len;
int num;
vector[num] returns[len];
corr_matrix[num] corr;
}
transformed data {
cholesky_factor_corr[num] L;
L = cholesky_decompose(corr);
}
parameters {
vector[num] mu;
vector<lower=-1, upper=1>[num] phi;
vector<lower=0>[num] sigma;
vector[num] h_std[len];
}
transformed parameters {
vector[num] h[len];
h[1] = mu;
for (t in 2:len) {
h[t] = h_std[t] .* sigma;
h[t] += mu;
h[t] += phi .* (h[t-1] - mu);
}
}
model {
mu ~ normal(0, 25);
phi ~ uniform(-1, 1);
sigma ~ inv_gamma(2.5, .025);
for (t in 1:len) {
h_std[t] ~ std_normal();
returns[t] ~ multi_normal_cholesky(rep_vector(0, num), diag_pre_multiply(exp(h[t] / 2), L));
}
}
```