Thanks again, Bob. Below is the code if any one else is interested.
generated quantities {
vector [L] y_new [T];
matrix[L,L] H_samp;
vector[L] mus_samp[T];
for (i in 1:T){
mus_samp[i] = (betafix[1] + alpha[i]) + X[i]*betafix[2];
}
for (i in 1:(L-1)) {
for (j in (i+1):L) {
H_samp[i,j] = sigma2 * exp(-d[i,j]*rho);
H_samp[j,i] = H_samp[i,j];
}
}
for (k in 1:L) {
H_samp[k, k] = sigma2 + tau2;
}
for (i in 1:T){
y_new[i] = multi_normal_rng(mus_samp[i], H_samp);
}
}