%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%***模型稳态求解
function [ys,check] =mythird_steadystate(ys,exo)
%%



global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

% initialize indicator
check = 0;
beta=0.99943;
alpha=0.55;  
mu  =	0.3; 
muzplus=muzplus_p;
muz=muz_p;
mupsi=1;
pistar=pitarget_p;
pitarget=pitarget_p;

y_obs=0;  
c_obs=0;  
i_obs =0;  
x_obs =0;  
m_obs =0;   
nu_obs =0; 
loan_obs=0;
r_obs =0; 
pic_obs =0;  
pii_obs=0; 
sbar_obs=1;
ystar_obs =0; 
pistar_obs =0;
rstar_obs=0;





pi=pitarget_p;
pic=pitarget_p;
pimx=pitarget_p;
pimc=pitarget_p;
pimi=pitarget_p;
pii=pitarget_p/mupsi;
piw=muzplus*pitarget_p;
pix=pistar;
s=1;
R=(muzplus_p/(beta)-taub)/(1-taub)*pitarget_p;
Rstar=(muzplus/beta-taub)/(1-taub)*pistar/(1+tau);
u=1;
Stilde=0;
Sprimetilde=0;
epsilon=1;
pitarget=pitarget_p;
Upsilon=1;
zetac=1;
phitilde=0;
zetah=1;

pitilded=pitarget_p^kappad*pitarget^(1-kappad);
mc=1/lambdad*(1-beta*xid*(pitilded/ pitarget_p)^(lambdad/(1-lambdad)))...
    /(1-beta*xid*(pitilded/ pitarget_p)^(1/(1-lambdad)))...
    *((1-xid*(pitilded/ pitarget_p)^(1/(1-lambdad)))/(1-xid))^(1-lambdad);
phalo=((1-xid*(pitilded/ pitarget_p)^(1/(1-lambdad)))/(1-xid))^(1-lambdad)...
    /((1-xid*(pitilded/ pitarget_p)^(lambdad/(1-lambdad)))/(1-xid))^((1-lambdad)/lambdad);

pitildex=pix^kappax*pix^(1-kappax);
mcx=1/lambdax*((1-xix*(pitildex/pix)^(1/(1-lambdax)))/(1-xix))^((1-lambdax))...
    *(1-beta*xix*(pitildex/pix)^(lambdax/(1-lambdax)))...
    /(1-beta*xix*(pitildex/pix)^(1/(1-lambdax)));
phalox=((1-xix*(pitildex/pix)^(1/(1-lambdax)))/(1-xix))^(1-lambdax)...
    /((1-xix*(pitildex/pix)^(lambdax/(1-lambdax)))/(1-xix))^((1-lambdax)/lambdax);

pitildemc=pimc^kappamc*pitarget^(1-kappamc);
mcmc=1/lambdamc*((1-ximc*(pitildemc/pimc)^(1/(1-lambdamc)))/(1-ximc))^((1-lambdamc))...
    *(1-beta*ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)))...
    /(1-beta*ximc*(pitildemc/pimc)^(1/(1-lambdamc)));
phalomc=((1-ximc*(pitildemc/pimc)^(1/(1-lambdamc)))/(1-ximc))^(1-lambdamc)...
    /((1-ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)))/(1-ximc))^((1-lambdamc)/lambdamc);

pitildemi=pimi^kappami*pitarget^(1-kappami);
mcmi=1/lambdami*((1-ximi*(pitildemi/pimi)^(1/(1-lambdami)))/(1-ximi))^((1-lambdami))...
    *(1-beta*ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)))...
    /(1-beta*ximi*(pitildemi/pimi)^(1/(1-lambdami)));
phalomi=((1-ximi*(pitildemi/pimi)^(1/(1-lambdami)))/(1-ximi))^(1-lambdami)...
    /((1-ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)))/(1-ximi))^((1-lambdami)/lambdami);

pitildemx=pimx^kappamx*pitarget^(1-kappamx);
mcmx=1/lambdamx*((1-ximx*(pitildemx/pimx)^(1/(1-lambdamx)))/(1-ximx))^((1-lambdamx))...
    *(1-beta*ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)))...
    /(1-beta*ximx*(pitildemx/pimx)^(1/(1-lambdamx)));
phalomx=((1-ximx*(pitildemx/pimx)^(1/(1-lambdamx)))/(1-ximx))^(1-lambdamx)...
    /((1-ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)))/(1-ximx))^((1-lambdamx)/lambdamx);

pitildew=pic^kappaw*pitarget^(1-kappaw) *muzplus;
whalo=((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw)...
    /((1-xiw*(pitildew/piw)^(lambdaw/(1-lambdaw)))/(1-xiw))^((1-lambdaw)/lambdaw);

pitildedpi=pitilded/pi;
pitildexpi=pitildex/pix;
pitildemcpi=pitildemc/pimc;
pitildemipi=pitildemi/pimi;
pitildemxpi=pitildemx/pimx;
pitildewpi=pitildew/piw;

%方程求解的稳态
    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xxx=[0.2856
     0.4940
     1.0267
    6.795  
    1.7784
    0.2500
    1.3551
    0.0398
    0.3241];
    
        opts=optimset('Display','off','MaxIter',200,'TolFun',1e-9,'TolX',1e-9);
        [yyy,fval]=fsolve(@mythird_steadystate_solvecore,xxx,opts,tauk,muzplus_p,alpha,lambdad,...
            etag,etax,omegax,etaf,omegai,etai,beta,b, sigmaL,omegac,...
            etac, etaa, lambdamc,taub,pitarget_p, delta,tau,...
          lambdami,lambdamx,lambdax,mu,Fomegabar,xiw,lambdaw,xid,xix,ximc,ximi,ximx,muz_p,...
            xy_target,nk_target,AL) 

        sigma=yyy(1);
        omegabar=yyy(2);
        Rk=yyy(3);
        kbar=yyy(4);
        wbar=yyy(5);
        H=yyy(6);
        ystar=yyy(7);
        rkbar=yyy(8);
        varphi=yyy(9);


%结合方程求解才能求解的稳态

pmx= Rstar*varphi/mcmx;
pmc= Rstar*varphi/mcmc;
pmi= Rstar*varphi/mcmi;
pc=(1-omegac+omegac*pmc^(1-etac))^(1/(1-etac));
pinvest=(1-omegai+omegai*pmi^(1-etai))^(1/(1-etai));
px=(R/(mcx*varphi))*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax));
pkprime=pinvest/Upsilon;
y=(phalo^(lambdad/(lambdad-1))*epsilon*(kbar/H/(mupsi*muzplus))^alpha)...
            /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H; 
nk=(pkprime*(1-Rk/R*((1-mu)*normcdf((log(omegabar)+sigma^2/2)/sigma-sigma)+omegabar*(1-normcdf((log(omegabar)+sigma^2/2)/sigma)))));
gamma=(1-0.001*y/nk/kbar)/((Rk-R-((normcdf((log(omegabar)+sigma^2/2)/sigma-sigma)+omegabar*(1-normcdf((log(omegabar)+sigma^2/2)/sigma)))-( (1-mu)*normcdf((log(omegabar)+sigma^2/2)/sigma-sigma)+omegabar*(1-normcdf((log(omegabar)+sigma^2/2)/sigma))))*Rk)*pkprime/nk+R)*(pitarget_p*muzplus);         
n=nk*kbar;
g=etag*y;
q=varphi/pc;
i=kbar*(1-(1-delta)/(muzplus*mupsi))/Upsilon;
x=px^(-etaf)*ystar;
c=1/(zetah*AL)*((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw*(1+sigmaL))...
            *(wbar*(1)/(lambdaw*(1)))*(1-beta*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)))...
            *(zetac/(muzplus-b)*(muzplus-beta*b)/(pc*(1)))/(1-beta*xiw*(pitildew/piw)^(1/(1-lambdaw)))*H^(-sigmaL);     
a=etaa*y;
k=kbar*u;
psizplus=(zetac*(muzplus-beta*b)/(muzplus-b))/(c*pc*(1));
im=omegai*(pinvest/pmi)^(etai)*i;
xm=omegax*(1/pmx*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax)))^etax...
            *phalox^(-lambdax/(lambdax-1))*px^(-etaf)*ystar;
cm=omegac*(pc/pmc)^(etac)*c;

Kd=lambdad*psizplus*y*mc/(1-beta*xid*(pitilded/pi)^(lambdad/(1-lambdad)));
Fd=psizplus*y/(1-beta*xid*(pitilded/pi)^(1/(1-lambdad)));
Kx=lambdax*psizplus*varphi*px*x*mcx/(1-beta*xix*(pitildex/pix)^(lambdax/(1-lambdax)));
Fx=psizplus*varphi*px*x/(1-beta*xix*(pitildex/pix)^(1/(1-lambdax)));
Kmc=lambdamc*psizplus*pmc*mcmc*cm/(1-beta*ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)));
Fmc=psizplus*pmc*cm/(1-beta*ximc*(pitildemc/pimc)^(1/(1-lambdamc)));
Kmi=lambdami*psizplus*pmi*mcmi*im/(1-beta*ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)));
Fmi=psizplus*pmi*im/(1-beta*ximi*(pitildemi/pimi)^(1/(1-lambdami)));
Kmx=lambdamx*psizplus*pmx*mcmx*xm/(1-beta*ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)));
Fmx=psizplus*pmx*xm/(1-beta*ximx*(pitildemx/pimx)^(1/(1-lambdamx)));
Kw=(zetah*H^(1+sigmaL))/(1-beta*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)));
Fw=(psizplus/lambdaw*H*(1)/(1))/(1-beta*xiw*(pitildew/piw)^(1+lambdaw/(1-lambdaw)));
smallh=whalo^(lambdaw/(1-lambdaw))*H;
aofu=0;


%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(ii,:));
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end