The Jacobian contains NaNs

Hi all,

I’ve set up the model to run in dynare. The steady state values are correct (no residuals). All parameter values (36 of them) are initialised. I tried to put different parameter values, but it didn’t change the outcome (and yes some of the parameter values are odd, because the priors are flat and the model will be estimated in later stages, but I am just trying to run it for now).

The error keeps croping up that “Jacobian contains NaNs or Inf”. I’ve tried using model_diagnostics command and it tells me “Input to SVD must not contain NaN or Inf”.

The steady state values of variables seem to make sense. I can’t think of any other angles so could someone help? As usual, the model contains non-standard features, but I’ve numbered the equations and described parameters and equations is some detail.

issue has been resolved

Equations 23 and 43 have terms like


The term in the bracket, evaluates to 0 given your steady state of w_HAT=W_HAT=1. Thus, a division by 0 occurs in the first derivative as it is raised to the power of phi_W_HAT-1-1=-.8

Thanks for your reply,

You are correct in that the term


evaluates to zero in the steady state (note: w_HAT=0, but W_HAT=1), but division by zero does not occur, because as you rightly calculated phi_W_HAT-1=0.2 (since the initialised value is equal to 1.2). And more importantly, even if phi_W_HAT was =1, the next term:


Multiplies the term that you’ve outlined. Both of them evaluate to 0 in the steady state and thus division by zero does not occur. So, unfortunately, I don’t think this is the problem.

However, I tried to use the initval command, and it seems to point to the steady state of hours:

Attempted to access params(38); index out of bounds because numel(params)=36.
Error in NL_ERPT_DSGE_steadystate2 (line 6)
    Error in evaluate_steady_state_file (line 54)
        [ys,params1,check] = h_steadystate(ys_init, exo_ss, params);
Error in evaluate_steady_state (line 58)
        [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M, ...
Error in resid (line 55)
    [oo_.steady_state,M_.params,info] = ...
Error in NL_ERPT_DSGE (line 666)
Error in dynare (line 180)
evalin('base',fname) ; 

As if there were 38 parameters that were entered in the steady state model block, though 36 of which were declared. I thoroughly checked both the steady state and the model blocks, but I can’t find any parameters that were not initialised. Do you have any further suggestions? Thank you.

Remember the derivative rules for exponents. When taking the derivative, decrease the exponent by 1. In the equation, you have phi_W_HAT-1=0.2 (since the initialised value is equal to 1.2). But in the Jacobian, you then get 0.2-1=-0.8. That was my typo above. For this reason, you get


which gives rise to the error.

Thank you. As you rightly suggested, the previous issue could be resolved simply by setting phi=>2, which turn out to be the only admissible values.

The code has slightly changed since the last time, but I have a new and seemingly less obvious issue. As usual, I don’t have any residuals, the steady state values seem to make sense. I am interested in solving this model using a second order approximation (or third) due to the non-linearities that are crucial in this paper. However, I find that I can only simulate the model using the first order approximation.

With higher order approximations, the correlations and autocorrelations can be computed, but the theoretical means are NaNs and without pruning IRFs are explosive (and not displayed). With pruning, IRFs are NaNs and not displayed. First order approximation seems to give rise to usual hump-shaped IRFs.

I am not quite sure what the problem is, but I would greatly appreciate some help.

issue has been resolved

It seems you shifted the problem exactly by one level. Now it is not the Jacobian, but the Hessian that contains the division by 0. Are you really sure that your equation is correct, i.e. that the base should evaluate to 0 in steady state?


The reason why I’m using this function is because I am trying to model downward wage rigidity (DELTA_W_HAT) and partially irreversible investment (DELTA_I), which are asymmetric adjustment cost functions. I could technically take out the term that is raised to the power PHI_I or PHI_W_HAT, but then the asymmetric component in the adjustment costs do not have a sufficient bite (its quantitatively small and makes little difference). Introducing this term was an idea of how to get an increase in the magnitude.

The reason I’m using this is because it gives rise to very neat first order conditions (as opposed to using some smooth transition function) and yet it can still be estimated using Bayesian methods (as opposed to using a piecewise-quadratic functional form).

In other words, there is a good reason why it is used, but I may have to come up with a slightly different functional form, unless you have any further suggestions?


Using a penalty function to embed a non-negativity constraint can be done. You might want to take a look at e.g. Kim/Kollmann/Kim (2010) Solving the incomplete market model with aggregate uncertainty using a perturbation method) or