# Need help: The steady state contains NaN or Inf

Dear all, I am new in Dynare, and I am trying to replicate the paper Backus, Kehoe, and Kydland (1994). “Dynamics of the trade balance and the terms of trade: the J-curve?”.

I wrote down 24 equations, 24 variables, the derivation steps of which can be seen from the document Appendix.pdf I attach to this post.

When I tried executing the codes in Dynare, Matlab keeps saying that “The steady state contains NaN or Inf”. I am really confused about what the problem is. I also attached the dynare codes to this post as well.

I really appreciate any kind of help! Thank you very much!

16. appendix 6.20.pdf (154.0 KB)

attempt.mod (2.1 KB)

This is a nonlinear model. Variables like output cannot be 0. So your initval-values do not make sense.

So am I supposed to change the initial values only? Or do I need to log-linearize the model first before entering them into Dynare?

Thank you very much!

No, if you linearize, you also need the steady states. Try to get the nonlinear model to work.

So do you mean I should manually calculate the steady states first?

If that is possible, that is always preferable. Use a nonlinear model with an analytically computed steady state

Thanks a lot for the reply!

I see. But a problem I have now is that when I tried calculating the steady states analytically, I got y* to be 0.

According to the stochastic process of shocks, I can write the law of motion of shocks as follows:

z1(+1) = 0.906 * z1 + 0.088 * z2 + \epsilon_1 (1)
z2(+1) = 0.088 * z1 + 0.906 * z2 + \epsilon_2 (2)

and the production function y1 = z1 * k1^(\theta) * n1^(1-\theta) (3).

Based on (1) & (2), if I set z1(+1) = z1, z2(+1) = z2, then the steady state of z1 and z2 are 0 (is that right?), which, based on (3), means the steady state of y1 is also 0. I know this result is weird, but I don’t know what the problem is.

Could you point out where I misunderstand things?

Thank you!

There must be a mistake in the original paper. It is standard to have mean 0 exogenous processes. But in that case either the production function must have the shock in exp()
y_t=\mathbf{e_z^z}k_{t-1}^\theta n_t^{1-\theta}
or the exogenous process is in logs

log(z1) = 0.906 *log(z1(-1)) + 0.088 *log(z2(-1)) + \epsilon_1
log(z2) = 0.088 *log(z1(-1)) + 0.906 *log(z2(-1)) + \epsilon_2


In both cases, the shock term in the production function is 1 in steady state. Note that in Dynare, you have to enter the two shock processes with a contemporaneous timing, not the t+1 in the original paper.

P.S.: It would be great if you could contribute your replication to https://github.com/JohannesPfeifer/DSGE_mod once you are done.

Thank you so much for the reply! Sure I definitely will once I get the solution.

Actually I am now stuck on solving for a system of equations to get the steady state. Can I ask by “analytically”, do you mean solving for the system using Matlab or some other software to solve for the roots to the system of equaionts? I am now left with eight equations and eight variables in order to solve for the steady state, and I tried using Matlab to solve it, but it never comes up with a result, since this process takes up a lot of memory. May I ask when you encounter this kind of complicated system of equations to solve for steady states, do you usually do it manually or using some kind of software?

Following are the eight equations I wrote down:
syms a1 a2 b1 b2 n1 n2 c1 c2
k1 = 1;
k2 = 1;
eqn1 = a1 + a2 - k1^0.36 * n1^0.64;
eqn2 = b1 + b2 - k2^0.36 * n2^0.64;

eqn3 = c1 + 0.025 * k1 - (0.85 * a1^(1/3) + 0.15 * b1^(1/3))^3;
eqn4 = c2 + 0.025 * k2 - (0.85 * b2^(1/3) + 0.15 * a2^(1/3))^3;

eqn5 = 1/0.99 - 1 + 0.025 - 0.36 * (n1/k1)^0.64 * (0.85a1^(1/3) + 0.15 * b1^(1/3))^2 * 0.85 * a1^(-2/3);
eqn6 = 1/0.99 - 1 + 0.025 - 0.36 * (n2/k2)^0.64 * (0.85
b2^(1/3) + 0.15 * a2^(1/3))^2 * 0.85 * b2^(-2/3);

eqn7 = 0.66c1 - 0.34 * 0.64 * (1-n1) * (k1/n1)^0.36 * (0.85a1^(1/3) + 0.15 * b1^(1/3))^2 * 0.85 * a1^(-2/3);
eqn8 = 0.66c2 - 0.34 * 0.64 * (1-n2) * (k2/n2)^0.36 * (0.85b2^(1/3) + 0.15 * a2^(1/3))^2 * 0.85 * b2^(-2/3);

[sol_a1, sol_a2, solb1, sol_b2, sol_n1, sol_n2, sol_c1, sol_c2] = solve(eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8)

You need to do more with pencil and paper. If you take a look at the calibration section, they impose a symmetric equilibrium, so you should be left with only half as many equations. Also, please do not plug in numbers for the parameters yet. It makes it impossible to see which equations are used.

Got you. And also, when I write “resid;” before steady, I found the following situation:

Residuals of the static equations:

Equation number 1 : 0.132
Equation number 2 : 0.132
Equation number 3 : NaN
Equation number 4 : NaN
Equation number 5 : 0
Equation number 6 : 0
Equation number 7 : 0
Equation number 8 : 0
Equation number 9 : 0
Equation number 10 : 0
Equation number 11 : -0.775
Equation number 12 : -0.775
Equation number 13 : 0
Equation number 14 : 0
Equation number 15 : 0.74234
Equation number 16 : 0.74234
Equation number 17 : 1.93
Equation number 18 : -12.2791
Equation number 19 : 0
Equation number 20 : -0.8
Equation number 21 : -0.2
Equation number 22 : 0
Equation number 23 : 0
Equation number 24 : 0

I am wondering whether all residuals are supposed to be 0 in a correct set-up? From the results above, there are some equations with nonzero residuals and two equations with NaN residuals, I know there must be something wrong with the equations that have NaN residuals, but for the nonzero residuals, can they happen in a correct RBC model?

If you analytical steady state is correct, then all residuals must be zero. But first focus on the NaN. They often indicate a fundamental problem like a variable forgotten to set or a missing exp().

Hi Professor Pfeifer,

Thanks a lot for your reply! After changing the initial values of n1 n2 from 1 to 0.5, all residuals become real numbers and there are no NaN residuals any more. But there are still nonzero residuals obviously as can be seen from below. And the error message I got is “Impossible to find the steady state. Either the model doesn’t have a steady state, there are an infinity of steady states, or the guess values are too far from the solution”. I think I should continue calculating the analytical steady states. But before doing that, I am still concerned that the equations I write are not exact. Are there any ways to check whether I have written all 24 equations correctly, without any missing equations or redundant equations? Now that I have derived all real-number residuals, does it mean that the equations I write guarantee a solution?

Residuals of the static equations:

Equation number 1 : -0.019719
Equation number 2 : 0.01911
Equation number 3 : -0.51129
Equation number 4 : -0.46319
Equation number 5 : 0.017945
Equation number 6 : 0.012897
Equation number 7 : 0
Equation number 8 : 0
Equation number 9 : 0
Equation number 10 : 0
Equation number 11 : 0
Equation number 12 : 0
Equation number 13 : 0
Equation number 14 : 0
Equation number 15 : 0.0035389
Equation number 16 : 0.063008
Equation number 17 : -0.033059
Equation number 18 : 1.2622
Equation number 19 : -0.18189
Equation number 20 : -0.23177
Equation number 21 : 0.1554
Equation number 22 : 0
Equation number 23 : 0
Equation number 24 : -0.59453

attempt.mod (2.6 KB)

Try the following initial values and then work from there:

var z1, z2, y1, y2, a1, a2, b1, b2, k1, k2, c1, c2, n1, n2, g1, g2, s1, s2, x1, x2, q1, q2, p, nx;

varexo epsilon1, epsilon2;

parameters beta, mu, gamma, theta, delta, J, sigma, rho, omega_im, w1, w2;
beta = 0.99;
mu = 0.34;
gamma = -1.0;
theta = 0.36;
delta = 0.025;
J = 1;
sigma = 1.5;
rho = -0.33333;
omega_im = 0.15;
w1 = 0.85;
w2 = 0.15;

model;
//intratemporal substitution between labor and consumption 1 &2
(1- mu)* c1 = mu * (1-n1)* (w1 * a1^(-rho) + w2 * b1^(-rho))^(-1/rho -1)* w1 * a1^(-rho - 1)*(1-theta)*y1/n1;
(1- mu)* c2 = mu * (1-n2)* (w1 * b2^(-rho) + w2 * a2^(-rho))^(-1/rho -1)* w1 * b2^(-rho - 1)*(1-theta)*y2/n2;

//Euler equation 3&4
c1^(mu * gamma -1) * (1-n1)^(gamma*(1-mu)) = beta * c1(+1)^(mu*gamma -1)*(1 - n1(+1))^(gamma *(1-mu))*((w1*a1(+1)^(-rho)+w2*b1(+1)^(-rho))^(-1/rho -1) * w1 * a1(+1)^(-rho -1) * theta * y1(+1)/k1(+1) + 1- delta);
c2^(mu * gamma -1) * (1-n2)^(gamma*(1-mu)) = beta * c2(+1)^(mu*gamma -1)*(1 - n2(+1))^(gamma *(1-mu))*((w1*b2(+1)^(-rho)+w2*a2(+1)^(-rho))^(-1/rho -1) * w1 * b2(+1)^(-rho -1) * theta * y2(+1)/k2(+1) + 1- delta);

//Production functions 5&6
y1 = z1 * k1^theta * n1^(1-theta);
y2 = z2 * k2^theta * n2^(1-theta);

//Evolution of tech shocks 7&8
log(z1) = 0.906 * log(z1(-1)) + 0.088 * log(z2(-1))+epsilon1;
log(z2) = 0.088 * log(z1(-1)) + 0.906 * log(z2(-1))+epsilon2;

//Evolution of government shocks 9&10
g1=0.95 * g1(-1);
g2=0.95 * g2(-1);

//Capital accumulation  11&12
k1 = (1 - delta) * k1(-1) + s1;
k2 = (1 - delta) * k2(-1) + s2;

//Market clearing conditions 13&14
a1 + a2 = y1;
b1 + b2 = y2;

//Goods market clearing 15&16
c1 + x1 + g1 = (w1 * a1^(-rho) + w2 * b1^(-rho))^(-1/rho);
c2 + x2 + g2 = (w2 * a2^(-rho) + w1 * b2^(-rho))^(-1/rho);

//Price determination 17&18
p = (w2/w1)* (a1/b1)^(1/sigma);
p = (w1/w2)* (a2/b2)^(1/sigma);
p = q2/q1;

//q1, q2 19 & 20
c1 + x1 + g1 = q1 * a1 + q2 * b1;
c2 + x2 + g2 = q1 * a2 + q2 * b2; // symmetric

//Investment 21 & 22
x1 = s1;
x2 = s2;

//net exports 23 & 24
nx = (a2 - p * b1)/y1;
end;

initval;
y1 = 1;
y2 = 1;
z1 = 1;
z2 = 1;
a1 = 0.2;
a2 = 0.8;
b1 = 0.8;
b2 = 0.2;
c1 = 0.2;
c2 = 0.2;
k1 = 1;
k2 = 1;
x1 = 0.8;
x2 = 0.8;
s1 = (1-delta)*k1;
s2 = (1-delta)*k2;
n1 = 0.62;
n2 = 0.62;
g1 =0;
g2 =0;
nx =-0.8;
p = 1;
q2 = 1;
q1 =1;
end;

shocks;
var epsilon1; stderr 0.00852;
var epsilon2; stderr 0.00852;
var epsilon1, epsilon2 = 0.258;
end;

resid;
check;

stoch_simul(order=1, hp_filter=1600);


Also note that

1. Your government spending shocks are currently not yet correct. They are mean 0, but should be mean 0.2 as per Table 3. You need something like
c1 + x1 + 0.2*exp(g1) = q1 * a1 + q2 * b1;

1. As is, you get a steady state, but you do not hit the calibration targets of Section III. \omega_1 and \omega_2 for example need to be set so that y=1. Only then will the government spending share be correct.

Thank you so much for your kind reply! I will look into these issues now. I really appreciate your help!

Hi professor Pfeifer,

After some trials and errors, I found only when w1 = 0.1 and w2 = 0.9 can I get a calibrated y1 to be 0.992253; alternatively, when I set w1 = 0.08, w2=0.92, the calibrated y1 becomes 1.01705. The values do not seem to be very realistic, but I tried almost all one-decimal numbers near 0.85 and 0.15 but the error message always shows that “impossible to find a steady state”. The previous two attempts seem to be the few that work. Do you think I can use “w1=0.1, and w2 = 0.9”?

For the government shocks, I am now assuming the government shocks to be 0, as is the case in their benchmark model.

Now I encounter the problem “The specified covariances for the structural errors are not consistent with the variances as they imply a correlation larger than ±1”. I don’t understand what problem this means. Look forward to your reply. Thanks a lot!

attempt2.mod (2.6 KB)

Hi professor Pfeifer,

I just realized that when I set w1, w2 to be the values I mentioned, the BK condition cannot be verified. Can I follow their attempt of w1=0.85, w2=0.15, which makes y1 = 0.937363 (is it too far from 1)?

Actually I am really confused by the error message “The specified covariances for the structural errors are not consistent with the variances as they imply a correlation larger than ±1”. What kind of errors does it imply? Look forward to your reply. Thank you very much!

Actually after deleting the command “var epsilon1, epsilon2 = 0.258”, I got the impulse response functions! Now I will try adding government shocks, and will keep you posted. Thank you so much!