but in each iteration, I need to get the derivatives both for the first and third order expansions. That is, for each parameter, I need to get two sets of derivatives (for 1st and 3rd order). Could you please let me know if/how I can do this.
Note: I have already tried to just change the order using options_.order = 3;, but I wasn’t successful.

Thank you very much jpfeifer, it works now. As you correctly guessed, I was using the wrong order. Just a quick follow up question, and please don’t bother if there is not a good answer for that. When I tried the following piece of code:

dynare afr0test
f_11_3 = [oo_.dr.ghx oo_.dr.ghu];
stoch_simul(var_list_);
f_11_test3 = [oo_.dr.ghx oo_.dr.ghu];
sum(vec(f_11_test3 -f_11_3)) %showing if the first order derivatives have changed

to see if running the same file from outside (i.e. using the stoch_simul) has changed the first-order derivatives, I got 2.14e-6. Admittedly, I have a lot of variables in the models (66 variables), which means each derivative is just changed a tiny bit, but do you know why there is a change in the first place?

When calculating the stochastic steady state as proposed by Johannes in another part of this forum (thanks Johannes!) I came across something which might be a bug? (or I’m overlooking something). I’m also using my .mod file for other stuff, therefore, I specifiy the approximation order in a parameter file. This worked fine up to order 2, however, no I’m encountrering the following error message:

The derivatives were not computed for the required order. Make sure that you used the right order option inside the stoch_simul command