function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)

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


%
load para_BenchmarkModel.mat;   
 
set_param_value('BETA',BETA);
set_param_value('MU',MU);
set_param_value('DELTA',DELTA);
set_param_value('PI_SS',PI_SS);
set_param_value('TAO_SS',TAO_SS);
set_param_value('G_SS',G_SS);
set_param_value('Z_SS',Z_SS);

load para_DECOMP.mat;

VARTHETA = para_DECOMP(1,1);   % The habit in consumption.
ETA      = para_DECOMP(1,2);   % Inverse Frisch Elasticity: Frisch Elasticity [0.25,0.53] see https://www.researchgate.net/publication/276137008_Estimates_of_the_Frisch_Elasticity_of_Labor_Supply_A_Review
CHI      = para_DECOMP(1,3);   % utility of money: 0.005 http://home.cerge-ei.cz/brazdik/download/miu.pdf
PSI      = para_DECOMP(1,4);   % disutility of labour
EPSILON  = para_DECOMP(1,5);   % Indexation
THETA    = para_DECOMP(1,6);   % Price elasticity, this is calibrated to target steady state mark-up 1/mc=1.5 (Leeper)
OMEGA    = para_DECOMP(1,7);   % Calvo probablity
%}

  g=log(G_SS);

  mc=log((THETA-1)*(1+PI_SS)/THETA*(1-OMEGA*BETA*(1+PI_SS)^(THETA*(1-EPSILON)))/(1-OMEGA*BETA*(1+PI_SS)^((THETA-1)*(1-EPSILON))));
  %mc=0;
  rr=1/BETA-1;
  r_k=1/BETA-1+DELTA;
  w=log((exp(mc)*exp(0)*MU^MU*(1-MU)^(1-MU)*r_k^(-MU))^(1/(1-MU)));
% closed-form expression for labor is not possible, so we need a helper function
  n0=1/3;
  n= BenchmarkModel_steadystate_helper(n0,VARTHETA,BETA,TAO_SS,w,PSI,DELTA,r_k, MU, g, ETA);
  k=log(exp(w)*exp(n)/r_k*MU/(1-MU));
  y=log(exp(0)*exp(n)^(1-MU)*exp(k)^MU);
  t=log(TAO_SS*exp(w)*exp(n));
  b=log((exp(t)-exp(g))/rr);
  lamda=log(-(1-VARTHETA)/(1-BETA*VARTHETA)/(1-VARTHETA)*(1-BETA*VARTHETA)/(exp(n)*(DELTA*exp(w)*MU/(r_k*(1-MU))-exp(0)*(exp(w)*MU/(r_k*(1-MU)))^MU)+exp(g)));
  c=log((1-VARTHETA)/(1-BETA*VARTHETA)/(1-VARTHETA)*(1-BETA*VARTHETA)/exp(lamda));
  pi=PI_SS;
  R=(1+pi)*(1+rr)-1;
  tao=TAO_SS;
  m=log((1+pi)/(1-BETA+pi)*CHI/exp(lamda));
  h=(pi/(1+pi)*exp(m));
  invest=log(exp(y)-exp(c)-exp(g));
  dy=0;
  dii=0;
  dg=0;
  db=0;
  dm=0;
  R_obs=0;
  pi_obs=0;
  tao_obs=0;
  y_hat=0;
  g_hat=0;
  m_hat=0;
  
                                                          
% save steady state values of variables
  MC_SS  = mc;
  Y_SS   = y; 
  B_SS   = b;
  R_SS   = R;
  M_SS   = m;

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