DSGE模型问题，DSGE模型收敛不了，稳态值求不出来，图也画不出来，麻烦大神帮我改进

...........................................................................DSGE模型程序
     我附上了模型的原文请大神帮我看看哪里出了问题，并且帮助我怎么改进模型最终求解出结果，画出脉冲响应图
     %模型中存在的问题1：总产出y是外生的，所有的方程左右两边y都约掉了；
     %模型中存在的问题2：价格p1和p2虽然是内生的，但无法求解具体表达式，p1和p2的稳态值我是通过matlab求解二元二次方程值拟合出来的；
var y y1 y2 p1 p2 L L1 L2 K K1 K2 tao taostar w w1 w2 r r1 r2 u u1 u2 z z1 z2 C I T fai H F aH aF ls pai A xi a1 a2;

varexo eA etaostar ea1 ea2 exi ey;
parameters omiga1 omiga2 sigama cita1 cita2  rouA rouxi roua1 roua2 routaostar rouy sigama1 sigama2  a b deta    ang beita gamma epseloH1 epseloH2 namdaH namdaF epseloF1 epseloF2 ustar ystar B;

parameters ys y1s y2s p1s p2s Ls L1s L2s Ks K1s K2s taos ws w1s w2s rs r1s r2s us u1s u2s zs z1s z2s Cs Is Ts  sigamaeA    sigamaetaostar  sigamaea1  sigamaea2  sigamaexi  sigamaey  fais Hs Fs aHs aFs lss taostars pais xis As a1s a2s;




omiga1=0.5;
omiga2=0.5;
sigama=0.96;  %   要拟合出p1和p2，    我这个模型中产品替代弹性在1附近才行，否则p1和p2会是复数
cita1=0.75;
cita2=0.25;
rouA=0.9;
rouxi=0.9;
rou1=0.9;
rou2=0.9;
routaostar=0.9;
rouy=0.9;
sigama1=0.45;
sigama2=0.55;
a=0.005;
b=0.14167;
deta=0.1;
ang=0.5;
beita=0.96;
gamma=0.55;
epseloH1=0.058;
epseloH2=0.083;
namdaH=1/0.015;
namdaF=1/3;
epseloF1=0.2;
epseloF2=0.06;
ustar=0.035;
ystar=120;
B=0.5;
sigamaeA =0.01;
sigamaetaostar=0.01;
sigamaea1=0.01;
sigamaea2=0.01;
sigamaexi=0.01;
sigamaey=1;
p1s=1;
p2s=1;
ys=ystar;
us=ustar;
taostars=0.2;
taos=taostars/(2-taostars);
rs=1/beita-1+deta;
Ks=(1-taos)*ys/(rs+((1-us)*(rs-deta)/(gamma*(1-us+us*ang)/(1-gamma)-ang*us-1+us)));
Cs=(rs-deta)/(gamma*(1-us+us*ang)/(1-gamma)-ang*us-1+us)*Ks*gamma/(1-gamma)*(1-us+u*ang);
Is=deta*Ks;
Ts=taostars*ys;
zs=1;
z1s=1;
z2s=1;
u1s=0;
ys=120;
y1s=(omiga1/p1s*(1-taos)/(1+taos))^(sigama)*ys;
y2s=(omiga2/p2s*(1-taos)/(1+taos))^(sigama)*ys;
As=1;
a1s=1/2;
a2s=1/2;
xis=1;
pais=0;
K1s=(1-cita1)^(sigama1/(1-sigama1))*(Ks/ys)^(sigama1/(sigama1-1))*1/(As*(1-a1s))*y1s*(1/(1-taostars))^(sigama1/(sigama1-1));
K2s=(1-cita2)^(sigama2/(1-sigama2))*(Ks/ys)^(sigama2/(sigama2-1))*1/(As*(1-a2s))*y2s*(1/(1-taostars))^(sigama2/(sigama2-1));
L1s=(cita1)^(sigama1/(1-sigama1))*((1-taostars)*ys-Ks/ys)^(sigama1/(sigama1-1))*1/(As*a1s)*y1s*(1/(1-taostars))^(sigama1/(sigama1-1));
u2s=((us*(cita2)^(sigama2/(1-sigama2))*((1-taostars)*ys-Ks/ys)^(sigama2/(sigama2-1))*1/(As*a2s)*y1s*(1/(1-taostars))^(sigama1/(sigama1-1)))+L1s)/((cita2)^(sigama2/(1-sigama2))*((1-taostars)*ys-Ks/ys)^(sigama2/(sigama2-1))*1/(As*a2s)*y1s*(1/(1-taostars))^(sigama1/(sigama1-1))+L1s);
L2s=(cita2)^(sigama2/(1-sigama2))*((1-taostars)*ys-Ks/ys)^(sigama2/(sigama2-1))*1/(As*a2s)*y1s*(1/(1-taostars))^(sigama1/(sigama1-1))/(1-u2s);
Ls=L1s+L2s;
ws=((1-taostars)*ys-Ks*rs)/(1-us)/Ls;
r1s=(1-cita1)*p1s*y1s^(1/sigama1)*((1-a1s)*As)^((sigama1-1)/sigama1)*(z1s*K1s)^(-1/sigama1);
r2s=(1-cita2)*p2s*y2s^(1/sigama2)*((1-a2s)*As)^((sigama2-1)/sigama2)*(z2s*K2s)^(-1/sigama2);
w1s=cita1*p1s*y1s^(1/sigama1)*(a1s*As)^((sigama1-1)/sigama1)*((1-u1s)*L1s)^(-1/sigama1);
w2s=cita2*p2s*y2s^(1/sigama2)*(a2s*As)^((sigama2-1)/sigama2)*((1-u2s)*L2s)^(-1/sigama2);
fais=ang*ws;
aHs=w1s*L1s/w2s/L2s/namdaH/(1+w1s*L1s/w2s/L2s/namdaH);
Hs=log(w1s)/log(w2s)*epseloH2/epseloH1;
aFs=(r1s/r2s)^(1/namdaF)/(1+(r1s/r2s)^(1/namdaF));
Fs=r1s/r2s*epseloF2/epseloF1;
lss=((1-taostars)*ys-Ks*rs)/(1-taostars)/ys;

model;
%1 总生产函数
exp(y)=(omiga1*(exp(y1))^((sigama-1)/sigama)+omiga2*(exp(y2))^((sigama-1)/sigama))^(sigama/(sigama-1));
%2 对第一部门边际求导
omiga1*(exp(y))^(1/sigama)*(exp(y1))^(-1/sigama)=exp(p1)*(1+exp(tao))/(1-exp(tao));
%3 对第二部门边际求导
omiga2*(exp(y))^(1/sigama)*(exp(2))^(-1/sigama)=exp(p2)*(1+exp(tao))/(1-exp(tao));
%4 第一部门生产函数
exp(y1)=exp(A)*((1-cita1)*((1-exp(a1))*(z1*K1(-1)))^((sigama1-1)/sigama1)+cita1*(exp(a1)*(1-u1)*L1)^((sigama1-1)/sigama1))^(sigama1/(sigama1-1));
%5 第二部门生产函数
exp(y2)=exp(A)*((1-cita2)*((1-exp(a2))*(z2*K2(-1)))^((sigama2-1)/sigama2)+cita2*(exp(a2)*(1-u2)*L2)^((sigama2-1)/sigama2))^(sigama2/(sigama2-1));
% 技术冲击
A=rouA*A(-1)+eA;
% 间接税率冲击
exp(taostar)=(1-routaostar)*taostars+routaostar*exp(taostar(-1))+etaostar;
% 第一部门技术进步偏向冲击
exp(a1)=(1-roua1)*0.5+roua1*exp(a1(-1))+ea1;
% 第二部门技术进步偏向冲击
exp(a2)=(1-roua2)*0.5+roua2*exp(a2(-1))+ea2;
%6 第一部门资本收益率
exp(r1)=(1-cita1)*exp(p1)*(exp(y1))^(1/sigama1)*((1-exp(a1))*exp(A))^((sigama1-1)/sigama1)*(exp(z1)*exp(K1(-1)))^(-1/sigama1);
%7 第二部门资本收益率
exp(r2)=(1-cita2)*exp(p2)*(exp(y2))^(1/sigama2)*((1-exp(a2))*exp(A))^((sigama2-1)/sigama2)*(exp(z2)*exp(K2(-1)))^(-1/sigama2);
%8 第一部门工资价格
exp(w1)=cita1*exp(p1)*(exp(y1))^(1/sigama1)*(exp(a1)*exp(A))^((sigama1-1)/sigama1)*((1-u1)*exp(L1))^(-1/sigama1);
%9 第二部门工资价格
exp(w2)=(cita2)*exp(p2)*(exp(y2))^(1/sigama2)*(exp(a2)*exp(A))^((sigama2-1)/sigama2)*((1-u2)*exp(L1))^(-1/sigama1);
%10 劳动总收入等于各部门劳动收入之和
exp(w)*(1-u)*exp(L)=exp(w1)*(1-u1)*exp(L1)+exp(w2)*(1-u2)*exp(L2);
%11 资本总收入等于各部门资本收入之和
exp(r)*exp(z)*exp(K(-1))=exp(r1)*exp(z1)*exp(K1(-1))+exp(r2)*exp(z2)*exp(K2(-1));
%12 劳动力市场出清
exp(L)=exp(L1)+exp(L2);
%13 资本市场出清
exp(K)=exp(K1)+exp(K2);
%14 第二部门的资本利用率为1
exp(z2)=1;
%15 第一部门的失业率为0
u1=0;
%16 总失业人数等于各部门失业人数之和
u*exp(L)=u1*exp(L1)+u2*exp(L2);
%17 总资本利用等于各部门资本利用之和
exp(z)*exp(K)=exp(z1)*exp(K1(-1))+exp(z2)*exp(K2(-1));
%18 税收等于税率乘以总产出
exp(T)=exp(taostar)*exp(y);
%19 体现了营业税在增值过程中重复征收的过程
exp(taostar)*exp(y)=exp(tao)*(exp(p1)*exp(y1)+exp(p2)*exp(y2)+exp(y));
%20 总产出等于总需求
exp(y)=exp(C)+exp(I)+exp(T)-exp(fai)*u*exp(L)+exp(K(-1))/exp(xi)*(a*exp(z)^2+(b-2*a)*exp(z)+a-b);
%21 总产出等于总收入
exp(y)=pai+exp(r)*exp(z)*exp(K(-1))+exp(w)*(1-u)*exp(L)+exp(T);
%22 利润pai全部用于补偿调整成本
pai=exp(K(-1))/exp(xi)*(a*exp(z)^2+(b-2*a)*exp(z)+a-b);
%23 资本积累方程
exp(K)=(1-deta)*exp(K(-1))+exp(xi)*exp(I);
%24 失业保障金是工资的一定比例
exp(fai)=ang*exp(w);
%  投资边际效率冲击
xi=rouxi*xi(-1)-exi;
%25 拉格朗日函数对L求导
(1-gamma)*exp(C)/gamma/exp(L)=exp(w)*(1-u)+exp(fai)*u;
%26 奥肯定律
(u-ustar)/ustar=B*(exp(y)-ystar)/ystar;
%27 拉格朗日函数对K求导
exp(C(+1))/exp(C)=beita*(exp(r(+1))*exp(z(+1))*exp(xi)-exp(xi)/exp(xi(+1))*(a*exp(z)^2+(b-2*a)*exp(z)+a-b)+exp(xi)/exp(xi(+1))*(1-deta));
%28 拉格朗日函数对z求导
exp(r)=(2*a*exp(z)+b-2*a)/exp(xi);
%29 城乡教育经费占比差距与工资差距的关系
exp(w1)/exp(w2)=namdaH*exp(aH)/(1-exp(aH))*exp(L2)/exp(L1);
%30 城乡人力资本水平差距与工资差距的关系
w1/w2=epseloH1/epseloH2*exp(H);
%31 城乡基础设施经费占比差距与资本收益率差距的关系
exp(r1)/exp(r2)=(exp(aF)/(1-exp(aF)))^(namdaF);
%32 城乡基础设施水平差距与资本收益率差距的关系
exp(r1)/exp(r2)=epseloF1/epseloF2*exp(F);
%33 劳动收入份额表达式
ls=exp(w)*(1-u)*exp(L)/(exp(w)*(1-u)*exp(L)+exp(r)*exp(z)*exp(K(-1)));
%34 总产出由外生冲击决定
y=(1-rouy)*ystar+rouy*y(-1)+ey;
end;
initval;
y=log(ys);
y1=log(y1s);
y2=log(y2s);
p1=log(p1s);
p2=log(p2s);
L=log(Ls);
L1=log(L1s);
L2=log(L2s);
K1=log(K1s);
K2=log(K2s);
tao=log(taos);
taostar=log(taostars);
w=log(ws);
w1=log(w1s);
w2=log(w2s);
r=log(rs);
r1=log(r1s);
r2=log(r2s);
u=us;
u1=u1s;
u2=u2s;
z=log(zs);
z1=log(z1s);
z2=log(z2s);
C=log(Cs);
I=log(Is);
T=log(Ts);
pai=pais;
fai=log(fais);
H=log(Hs);
F=log(Fs);
aH=log(aHs);
aF=log(aFs);
ls=lss;
A=log(A);
xi=log(xis);
a1=log(a1s);
a2=log(a2s);
end;

resid;
steady;

shocks;
var  eA = sigamaeA^2;      
var etaostar = sigamaetaostar^2;              
var ea1 = sigamaea1 ^2;
var  ea2 = sigamaea2^2;
var exi = sigamaexi^2;
var  ey= sigamaey^2;
end;
stoch_simul(order=1);...........................................................................................回归结果非常不好



Residuals of the static equations:

Equation number 1 : 44.1053 : 1
Equation number 2 : 0 : 2
Equation number 3 : 130.6226 : 3
Equation number 4 : 75.8947 : 4
Equation number 5 : 75.8947 : 5
Equation number 6 : NaN : A
Equation number 7 : 0 : 7
Equation number 8 : NaN : 8
Equation number 9 : NaN : 9
Equation number 10 : -Inf : 10
Equation number 11 : -Inf : 11
Equation number 12 : -Inf : 12
Equation number 13 : -Inf : 13
Equation number 14 : -17018.2435 : 14
Equation number 15 : -517.4791 : 15
Equation number 16 : 0 : 16
Equation number 17 : -40.7392 : 17
Equation number 18 : 0 : 18
Equation number 19 : 0 : u1
Equation number 20 : -2.3585 : 20
Equation number 21 : -40.7392 : 21
Equation number 22 : 0 : 22
Equation number 23 : -6.1988 : 23
Equation number 24 : -9.5238 : 24
Equation number 25 : 46.2361 : 25
Equation number 26 : 0 : pai
Equation number 27 : -32.6372 : 27
Equation number 28 : 0 : 28
Equation number 29 : 0 : xi
Equation number 30 : 1.926 : 30
Equation number 31 : 0 : 31
Equation number 32 : 0 : 32
Equation number 33 : -3.3333e-06 : 33
Equation number 34 : 0 : 34
Equation number 35 : 0 : 35
Equation number 36 : 0 : 36
Equation number 37 : 0 : 37
Equation number 38 : -0.48025 : ls
Equation number 39 : -11.5213 : y


警告: Some of the parameters have no value (roua1, roua2) when using steady. If these parameters are not initialized in a steadystate file or a
steady_state_model-block, Dynare may not be able to solve the model. Note that simul, perfect_foresight_setup, and perfect_foresight_solver do not
automatically call the steady state file.
Randomize initial guess...




Residuals of the static equations:

Equation number 1 : NaN : 1
Equation number 2 : NaN : 2
Equation number 3 : NaN : 3
Equation number 4 : NaN : 4
Equation number 5 : NaN : 5
Equation number 6 : NaN : A
Equation number 7 : NaN : 7
Equation number 8 : NaN : 8
Equation number 9 : NaN : 9
Equation number 10 : NaN : 10
Equation number 11 : NaN : 11
Equation number 12 : NaN : 12
Equation number 13 : NaN : 13
Equation number 14 : NaN : 14
Equation number 15 : NaN : 15
Equation number 16 : NaN : 16
Equation number 17 : NaN : 17
Equation number 18 : NaN : 18
Equation number 19 : NaN : u1
Equation number 20 : NaN : 20
Equation number 21 : NaN : 21
Equation number 22 : NaN : 22
Equation number 23 : NaN : 23
Equation number 24 : NaN : 24
Equation number 25 : NaN : 25
Equation number 26 : NaN : pai
Equation number 27 : NaN : 27
Equation number 28 : NaN : 28
Equation number 29 : NaN : xi
Equation number 30 : NaN : 30
Equation number 31 : NaN : 31
Equation number 32 : NaN : 32
Equation number 33 : NaN : 33
Equation number 34 : NaN : 34
Equation number 35 : NaN : 35
Equation number 36 : NaN : 36
Equation number 37 : NaN : 37
Equation number 38 : NaN : ls
Equation number 39 : NaN : y


错误使用 print_info (line 32)
The steady state has NaNs or Inf.

出错 steady (line 102)
    print_info(info,options_.noprint, options_);

出错 lsdsge.driver (line 777)
steady;

出错 dynare (line 278)
    evalin('base',[fname '.driver']);

Caught MathWorks::System::FatalException