Correlated errors in a factor model

Consider a standard factor model such as the one taken from this blog:

data {
  int<lower=1> n;
  int<lower=1> nX;
  vector [n] y;
  matrix [n,nX] X;

parameters {
  vector[nX] beta;
  real<lower=0> sigma;

transformed parameters {
  vector[n] mu;
  mu = X*beta;

model {
  // Likelihood
  // Priors
  beta ~ normal(0,100);

generated quantities {
  vector[n] log_lik;
   for (i in 1:n) {
   log_lik[i] = normal_lpdf(y[i] | mu[i], sigma); 

Now consider that one wants to add correlated residuals (error terms) between items in this model. What would be the best way to do so? Do you need to literally calculate the residuals in the generated quantities block and then assume a multivariate normal distribution for them in the model block? Or should this be specified in some other way?

John, I don’t think it is possible to model things that are defined in generated quantities. If you are looking to have just a few correlated residuals, you might be able to add something to mu as shown here (scroll down to the “Incorporating lagged residuals” section).

If you want all the residuals to be correlated, you probably have to move to a multivariate normal distribution, where sigma goes down the diagonal of the covariance matrix and off-diagonal entries are for your residual covariances.