Could you please expand a bit on why things could get tricky? Or if you could provide me with any reference, I shall be very happy to check them out myself!

I am not looking for the optimal policy function, but nonetheless I am having trouble solving the System as you have done in the cited code, with the opportune modifications.

In particular, I am now running the latest version of Dynare in which the Homotopy approach works. The result of such method is “Failed to solve for the perfect foresight model”, i.e. it does not find convergence (Although it produces a simulation which would go accordingly with my theory).

I have tried by reducing the size of the shock, but I think I do need a bigger one for it to have the desired effects. Therefore, I was wondering if anyone could suggest an alternative strategy (or again, provide a useful reference).

For smaller shocks (up to a drop in real interest rate to -0.30 from the initial value of 1) the Homotopy Method does work, but there is no discernible difference with respect to the standard, quadratic model.

If I set a drop of the real interest rate to -1, however, the homotopy method continues for around 1300 steps with progressively increasing lambdas, but also increasing max residuals (in the order of e+04, which is really huge). Some of them succeed, some do not. In particular the first step (Lambda=0.5) has very small residuals, which makes it try for Lambda=1, from which the divergence starts.

I am not quite sure on how this residuals are being calculated, though. Maybe understanding how these residuals are being calculated is a first step towards understanding the nature of my issue.

Each equation at each point in time has a residual. What is displayed is the infinity norm (maximum absolute value) of these residuals. Are you sure that the nonlinear model has a solution with big shocks?

In fact, I am not sure whether the model does have a solution for big shocks (it might be argued that a natural interest rate of -1 is not such a big shock, being the standard one used by Galì). That is what I was hoping to figure out.

That seems to be a good place to start. Do you mean setting the initial values of the simulation equal to the values of the solution after the shock?

Otherwise I was thinking about moving to another method, always within dynare, the most famous of which is probably OccBin by Iacoviello. However, I believe that it works only for linear models, as most other methods do.

Thank you so much for any help you can provide me!

I will do that straight away and hopefully that will provide me some further answers.

I was thinking that maybe it would be appropriate to generate multiple simulations for different values of the key parameter and the key shock, showing where the model does get solved and where not. Is there any particular reference on how to implement such an analysis, by inserting a Dynare routine within a more general .m file?

Thank you so much again for the quick replies and the help!

I have a couple more query. I am sorry to ask you so many things, but you have been very helpful!

I was also thinking of trying to apply the extended path method to the model, to have an approximate result of what would happen if expectations were to be taken into account.

Even without the ZLB, I run into an error message when I try to simulate an exogenous cost-push shock, as follows

[code]
%define natural rate shock
shocks;
% var r_nat;
% periods 1:6;
% values -1;
var u_cp; stderr 0.01;
end;

extended_path(periods=100,order=2);[/code]

In particular, the error is

Is there any way in which you could provide me a reference where I could see how the extended path method is to be implemented in Dynare, e.g. some explanatory mod files? If they included the ZLB, then it would be even more awesome.

Do you think it would be possible to simulate the model stochastically with a cost-push shock via stoch_simul(periods=300, order=2), even though the model presents a nonlinear target function, in particular x_t=-k*(exp^(a*pi)-1)/a ? Clearly, without ZLB. I have done it and the model gets solved and changes for different parametrization of a, but I wanted to understand whether the results make sense.

As long as your nonlinear function is differentiable, that should work. Even if it is not differentiable, you can still go for a penalty function approach. See the last post at [The Jacobian contains NaNs)

I have been trying to apply a penalty function approach lately, so as to implement the ZLB without the need of the max function in a stochastic context, but dynare seems to have some difficulty in finding the steady state. I have tried to lower the value of the variance of the shock, but I had no luck there.

I have the same issue but unfortunately the code you posted in your response is no longer visible. Would you mind to once again show how to use the solution from the linear version as a starting point the nonlinear estimation.