Non Centered Parameterization in GP Model

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);
  }
}