function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% 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
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright (C) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

betta =0.998;
phi=1.17;
gamma=0.7;
psi=1;
xi =0.82;
lambdaf=10/9;

%% Enter model equations here

options=optimset();   % set options for numerical solver

PI=PIbar;             % 通胀内生变量 = 稳态的通胀参数值
z =1;
R=PI/betta;           % 名义利率，由于内生变量 PI 和 mu_z 都已经表示成了参数的函数，此处相当于 R 也表示成了参数的函数
ptilt = ((1-xi*(PI)^(1/(lambdaf-1)))/(1-xi))^(1-lambdaf); 
pstar = (1-xi*(PI)^(lambdaf/(lambdaf-1)))/((1-xi)*ptilt^(lambdaf/(1-lambdaf)));
s= (ptilt/lambdaf)*(1-betta*xi*PI^(lambdaf/(lambdaf-1)))/(1-betta*xi*PI^(1/(lambdaf-1)));
w=(s/((1/(1-gamma))^(1-gamma)*(1/gamma)^gamma*(1-psi+psi*R)))^(1/gamma);

tempHtoI = (1/w)*(gamma/(1-gamma));
tempYtoI = pstar*tempHtoI^gamma;
tempCtoY = 1- (1/tempYtoI);

h = (gamma/(tempYtoI*tempCtoY*(1-gamma)))^(1/(1+phi));
ivst = h/tempHtoI;
y = tempYtoI*ivst;
c= tempCtoY*y;
lambda= 1/c;
X1 = (lambdaf*lambda*y*s)/(1-betta*xi*PI^(lambdaf/(lambdaf-1)));
X2 = (lambda*y)/(1-betta*xi*PI^(1/(lambdaf-1)));

%% end own model equations

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

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