function [ys,params,check] = benchmark_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)

% 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

% Technology
a = a_ss;

% Tax
tau= tau_ss;

% Chi
chi = chi_ss;

% Gamma 
gamma = gamma_ss; 

% Eqn 2: Bond Euler
r  = 1 / beta_p;

% Eqn 7: Investment Demand
b_k = ( (eta_p * tau * r) / (nu_p * omega_p) )^( 1/ (omega_p - 1) ); 


% Eqn 8: Capital Demand
aux1 =(1 - tau * chi - gamma) * (1/eta_p - 1 + delta_p);
aux2 = eta_p * tau * r * b_k ;
aux3 = nu_p * b_k^omega_p; 

y_k = (aux1 - aux2  + aux3) * 1/ ((1-tau)) * 1/alpha_p;

% Eqn 4: Production Function
h_k = (y_k /a )^ (1 / (1 - alpha_p ) );


% Eqn 6: Labor Demand
w = (1 - alpha_p) * a * h_k^(-alpha_p);

% Eqn 5: Capital accumulation
i_k = delta_p;

% Eqn 3: Resource constraint
c_k = y_k - i_k; 

% Eqn 1: Consumption labor choice
k = ( w / theta_p * h_k ^ (-phi_p) * c_k ^ ( -sigma_p )) ^ (1/( sigma_p + phi_p ));

y = y_k * k;
c = c_k * k;
i = i_k * k;
h = h_k * k;
b = b_k * k;

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