If you set order=2 and simulate a long enough time series, Dynare will provide you with E. Given that you know the functional form for utility as well as the steady state values for c and h, you can analytically solve
\frac{((1-\xi)*c_ss)^{1-\sigma}}{1-\sigma}-\varphi\frac{n_ss^{1+\nu}}{1+\nu}=E
for \xi.
You might also be able to use theoretical moments, but there one would need to check whether second moments based on the first order approximation a la Kim/Kim/Schaumburg/Sims are sufficient.