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

% 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

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

[(ce)^(-psi)]=(1+((q)-1)/(1-etat*(q)))*((cw)^(-psi))*(1-(N))^(omega*(1-psi));
[omega*(cw)/(1-(N))]=w;
q=beta*(demt(+1)/demt)*gss*((cw(+1))/(cw))^(-psi)*((1-(N(+1)))/(1-(N)))^(omega*(1-psi))*(R(+1)*(1))*(q(+1))*(1-delta)+theta*(((q(+1))-1)/(1-eta(+1)))*(q(+1))*((R(+1))*(1))+(pi(+1))*(q(+1))*(1-delta);
[1]=beta*(demt(+1)/demt)*((cw(+1))/(cw))^(-psi)*((1-(N(+1)))/(1-(N)))^(omega*(1-psi))*(i);
[ce+(1-etat*(q))*(I)]=((R)*(1)*(k(-1)))+(zeta)*(q)*((1-delta)*(k(-1)));
y=(exp(Z))*(k(-1))^alpha*((1-zeta)*(N))^(1-alpha)- phi;
y=theta*(ce)+(1-theta)*(cw)+zeta*(I)+G+(((mu)/(2))*(pi- 1)^2*(y));
k=(1-delta)*(k(-1))+zeta*(I);
[k/N]=((alpha)/(1-(alpha)))*((w)/(R));
[(mu)*((pi)-1)*((pi))]=(mu)*beta*(lamda(+1))*((pi(+1))-1)*(pi(+1))*((y(+1))/(y))+1-v+(v)*MC+((mu)/(2))*(v)*(((pi(+1))-1)^2);
MC=(1-alpha)^-(1-alpha)*(alpha)^-(alpha)*(w)^(1-alpha)*(R)^alpha*(1/Z);
cspt=theta*(ce)+(1-theta)*(cw);
zeta=sszeta;
gammat=((q)-1)/(1-(q)*etat);

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