Hi @Elef.
A couple of points: your posted Stan code does not include a squared term for x (i.e. you need beta_sq * square(x[n]) in the model formula) and your parameters are varying by data case, whereas the model you are trying to simulate has, I believe, constant parameters (i.e. you want normal_rng(alpha + beta * x[n] + beta_sq * square(x[n]), sigma) as the linear predictor.