各位好，
我现在在尝试replicate Backus, Kehoe, and Kydland 1994年的paper：Dynamics of the trade balance and the terms of trade: the J-curve? 我在下面附上了我写的Dynare codes，24个式子，24个变量。持续弹出的问题如下：
Error using print_info (line 90)
The steady state contains NaN or Inf
Error in steady (line 104)
print_info(info,options_.noprint, options_);
Error in attempt (line 266)
steady;
Error in dynare (line 223)
evalin('base',fname) ;


我的代码如下：

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- 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

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);

y1 = z1 * k1^theta * n1^(1-theta);

y2 = z2 * k2^theta * n2^(1-theta);

z1 = 0.906 * z1(-1) + 0.088 * z2(-1)+epsilon1;

z2 = 0.088 * z1(-1) + 0.906 * z2(-1)+epsilon2;

g1=0.95 * g1(-1);

g2=0.95 * g2(-1);

k1 = (1 - delta) * k1(-1) + s1;

k2 = (1 - delta) * k2(-1) + s2;

a1 + a2 = y1;

b1 + b2 = y2;

c1 + x1 + g1 = (w1 * a1^(-rho) + w2 * b1^(-rho))^(-1/rho);

c2 + x2 + g2 = (w2 * a2^(-rho) + w1 * b2^(-rho))^(-1/rho);

p = (w2/w1)* (a1/b1)^(1/sigma);

p = (w1/w2)* (a2/b2)^(1/sigma);

p = q2/q1;

c1 + x1 + g1 = q1 * a1 + q2 * b1;

c2 + x2 + g2 = q1 * a2 + q2 * b2; // symmetric

x1 = s1;

x2 = s2;

nx = (a2 - p * b1)/y1;

end;

initval;

y1 = 1;

y2 = 1;

z1 = 0.1;

z2 = 0.1;

a1 = 0.1;

a2 = 0.9;

b1 = 0.1;

b2 = 0.9;

c1 = 0.8;

c2 = 0.8;

k1 = 10;

k2 = 10;

g1 =0;

g2 =0;

end;

shocks;

var epsilon1; stderr 0.00852;

var epsilon2; stderr 0.00852;

var epsilon1, epsilon2 = 0.258;

end;

steady;

check;

stoch_simul(order=1, hp_filter=1600);

好像是steady states里全部是NaN和Inf。题目里的信息是这样的：
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)

在计算steady states的时候对于techonology shock zt我是这么求的：让z1(+1) = z1, z2(+1) = z2, 这样的话求出来的z1 z2 steady state上都是0，那么y1* y2*也都是0了。感觉不大对。另外需要手动把24个方程log-linearize吗？实在是非常需要各位的建议！感谢啊！

我在帖子后也附加了求euler equation和intratemporal substitution那几个式子的过程，各位如果有时间也可以过目一下。感激不尽！