function [ys,params,check] = RBC_MPShock_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = RBCSIMS_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)



alpha=0.33;
beta=0.99;

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


%% Enter model equations here

%mu_ss = 0;

q = 1;
ttau_ss = 0.2;
P_ss = varepsilon*(1-alpha)/(1-varepsilon);
mu_ss = (((1-alpha)/psiss)*((1/beta)- (1-delta))-alpha*(1-ttau_ss));
kn_ss = (PP*alpha*(1-ttau_ss)/(1/beta -(1-delta)-mu_ss*psiss))^(1/(1-alpha));
n = ((1/theta)*((1-ttau_ss)/(1 - ttau_ss + mu_ss))*(1/(kn_ss^alpha - delta*kn_ss))*(1-alpha)*(kn_ss^alpha))^(1/(1+chi));
k = kn_ss * n;
psi = psiss;
nint = (1/beta) - 1;
y = k^alpha*n^(1-alpha);
i = delta*k;
c = y - i;
lambda = 1/(c*PP);
w = psiss*kn_ss;
PP = P_ss;
mu = mu_ss;
ttau = ttau_ss;


A = 1;
log_c=log(c);
log_y=log(y);
log_k=log(k);
log_n=log(n);
log_i=log(i);
log_w =log(w);
log_r = log(r);
log_A = log(A);
log_psi = log(psi);
log_P = log(PP);


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