This is likely the same issue you encountered in your other thread here, because you have a hierarchical model.
That is, the hyperparameters for theta (mu and tau) are themselves parameters:
target += normal_lpdf(theta[j,k] | mu[j],tau[j]);
Have a look at the User’s guide chapter that I linked in that topic and its discussion of the non-centered parameterisation for these models.
Also, if thetanew and ynew are predictions from the model, rather than model parameters, then it’s more efficient to compute those in the generated quantities block