Problem in replicating results with moments_varendo?

You can use a posterior_function. See Posterior predictive moments