Dear friends, sorry to post a simple problem. When I compute IRF at third order using Dynare routine, the simulated GIRF seems unstable. They change with such options as "periods = " , " irf = ". 
Will this problem be solved if I follow such as Prof. Jesus, Pfeifer’s method about computing IRF ? Thanks in advance!
If you use standard IRFs in Dynare at order=3, you should use many replications and pruning. Computing them at the stochastic steady state/ergodic mean in the absence of shocks (EMAS) is more stable, because it does not require and stochastic simulations. GIRFs with the Andreasen et al. toolkit are also more stable, because they are computed without simulations as well.
Dear professor, so much thanks for your reply and useful advice!
(Could I ask you another question? After approximation, what’s the meaning of results stored in oo .dr.g0, oo .dr.g1, oo .dr.g2, and oo .dr.g3? According to the manual, they should mean different kinds of the Kronecker product of state variables except for oo .dr.g0.)
Yes, those matrices store the coefficients of the decision rules associated with the steady state, first, second, and third order terms.
Hi, professor! Thanks for your reply and help. But I am still confused about why the second-moments are usually large while the IRFs computed by EMAS are so small? How can I examine the effects of volatility shocks on moments? In the appendix of Pro. Pfeifer’s paper (2014), it exhibits first-moments of variables in difference. 
Thanks for your help in advance!
What do you mean? Most of the movements in the data come from the level shocks, not the volatility shocks.
Hi, professor, thanks for your reply and advice!
But why the moments computed at higher order approximation are usually larger than that at first order? Thanks a lot!