I did two ways. But appears both are wrong.

set z1=x(+1)y(+1), z2=x(+1), and use z1/z2

x(+1)y(+1)/z2
I compared the IRF of 1, 2 with the result if I use x(+1)y(+1)/x(+1). They ended the same. However, I think x(+1)y(+1)/x(+1) is treated as E x(+1)y(+1)/x(+1)] in Dynare. I am pretty sure there is some correlation between x and y, so that E(x(+1)y(+1))should not equal E[x(+1)]E[y(+1)]. I was expecting different results between 1 or 2 and x(+1)y(+1)/x(+1).
What is the correct way to express E(x(+1)y(+1))/E(x(+1)) in Dynare? Cheers!
At first order, there will be no differences at all. You need to go at least to second order to see any differences.
Both ways, 1 and 2 with the auxiliary variables should be the same, although I would prefer the first one.
Thank you. Can you explain why the difference shows up with the orders above 2?
Actually, it turned out when I put in stoch_simul(order = 3). How is it decided? Maybe the same question.
The reason is that first order approximations result in certainty equivalence. Due to linearization only first order terms appear. But covariances are second order terms and thus only show up at second order or higher.