%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%***模型稳态求解
function [ys,check] =mysecond_steadystate(ys,exo)
% function [ys,check] = myfirst_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)


global M_ options_

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


%% Enter model equations here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%***直接确定的稳态，等于1或者某个固定的参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%***外生冲击变量，为1或者参数值，8个,sigma需要计算,g后面给出
zstar=1;
zetac=1;
zetah=1;
pitarget=pitarget_p;
muz=muz_p;
up=1;
varepsilon=1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%***核心变量中直接确定的稳态，17+2个
s=1;
wbar=1;
pbar=1;
pxbar=1;
pmcbar=1;
pmibar=1;
pmxbar=1;

pi=pitarget_p;
pic=pitarget_p;
pii=pitarget_p;
pix=pitarget_p;
pimc=pitarget_p;
pimi=pitarget_p;
pimx=pitarget_p;
pistar=pitarget_p;

nfa=0;

u=1;

rd=pitarget_p*muz_p/beta_p;



rstar=pitarget_p*muz_p/(beta_p*(1+taub_p));

sbar=1;
gammac=(1-omegac_p+omegac_p*((1-omegac_p)/((sbar*lambdamc_p)^(etac_p-1)-omegac_p)))^(1/(1-etac_p));
gammai=(1-omegai_p+omegai_p*(lambdami_p*sbar*gammac)^(1-etai_p))^(1/(1-etai_p));
gammax=lambdax_p*(1-omegax_p+omegax_p*(lambdamx_p*sbar*gammac)^(1-etai_p))^(1/(1-etai_p))/(sbar*gammac);

rk=(rd/pitarget_p-1+delta_p)*gammai  ;
w=rk*(1-alpha_p)/(alpha_p*muz_p)*1/((lambdad_p*rk*muz_p^(alpha_p-1)/alpha_p)^(1/(1-alpha_p)));
%%%%%%
%稳态核心计算的实验
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%初值
x0=[1
    0.3
    0.7
    0.3
    1.3];
%%
opts=optimset('Display','off','MaxIter',2000,'TolFun',1e-20,'TolX',1e-20);


[yyy,fval] = fsolve(@mysecond_steadystate_core,x0,opts,omegac_p,etac_p,lambdamc_p,omegai_p,etai_p,...
                lambdami_p,omegax_p,etax_p,lambdamx_p,lambdax_p,lambdad_p,muz_p,alpha_p,...
            delta_p,etaf_p,al_p,sigmal_p,lambdaw_p,b_p,beta_p,pitarget_p)

k=yyy(1);
h=yyy(2);
y=yyy(3);
c=yyy(4);
ystar=yyy(5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%***依据上述核心求解计算的，13个
psiz=(1/((1-b_p/muz_p)*c)-beta_p*b_p/((muz_p-b_p)*c))/gammac;

q=gammai;

kbar=k;

i=(1-(1-delta_p)/muz_p)*k;

gammamc=lambdamc_p*sbar*gammac;
gammami=lambdami_p*sbar*gammac;
gammamx=lambdamx_p*sbar*gammac;

fd=psiz*y/(1-beta_p*xid_p);
fw=psiz*h/(lambdaw_p*(1-beta_p*xiw_p));
fmc=psiz*gammamc*omegac_p*(lambdamc_p*sbar)^(-etac_p)*c/(1-beta_p*ximc_p);
fmi=psiz*gammami*omegai_p*(lambdami_p*sbar)^(-etai_p)*i/(1-beta_p*ximi_p);
fmx=psiz*gammamx*omegax_p*((omegax_p*gammamx^(1-etax_p)+1-omegax_p)^(1/(1-etax_p))/gammamx)^etax_p*gammax^(-etaf_p)*ystar/(1-beta_p*ximx_p);
fx=psiz*sbar*gammac*gammax^(1-etax_p)*ystar/(1-beta_p*xix_p);

g=0.2 *y;    




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