I found where the problem originated.
At about line 84 of th_autocovariances.m you find the block

[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
iky = iv(ivar);
if ~isempty(u)
iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
ivar = dr.order_var(iky);
end

What happens is that lyapunov_symm returns a ‘u’ matrix that in the position of the selected variable (i.e. the variable(s) listed after stoch_simul) has some values that violate the Schur_vec_tol. It then returns an empty iky vector, which will produce an ivar empty vector etc.

Notice though that if I don’t list any variable after stoch_simul, Dynare returns a solution. with numeric (second order) constant term and numeric correction, although the mean vector is NaN.
Furthermore, as mentioned in the previous message in this thread, I can compute the (conditional and unconditional) mean of the variables to second order using the state-space solution that Dynare provides even in the buggy case.

The source of the problem seem to be the existence of a unit root in some (not all) variables.

I think that if this is correct a warning message in the block quoted above would be very much appreciated