I want to code in Dynare an equation that can be simplified to something as:
aE[X(t+1)]=bE[X(t+1)*Y(t+1)]
where E[X(t+1)*Y(t+1)] is the expectation of a product of two random variables

Any hint on how can I code this in Dynare?

if I define the variable Z=X*Y and I write:
aX(+1)=bZ(+1)
then Dynare seems to understand:
aE(X(t+1))=bE[X(t+1)]*E[Y(t+1)] , that is, a=bE[Y(t+1)]
but this is not what I want to write.

If you use order=1, this is exactly what you should get. In a linearized version, certainty equivalence holds and the covariances drop out. The trick with the auxiliary variable is usually the correct way for higher order approximations. If it does not work, please post the mod-file.

in previous posts I read Dynare automatically adds an E_t in front of all equations.
Thus, if I define the variable Z=X*Y, when I write z(+1), then Dynare understands E[X(t+1)*Y(t+1)], not E[X(t+1)]*E[Y(t+1)]. Is this correct?
And a different issue is that up to the 1st order the covariances/variances terms disappear