function [ys,check] = model_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% 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
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(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


Lambda-beta;
    
M - (m_c^((eta - 1)/eta) + theta*m_d^((eta - 1)/eta))^(eta/(eta - 1));
   
(Lambda - 1)*p + Lambda*D + psi_e*e^(-gamma_e)*pi^(1 - gamma_e);

c_e - D*e;

1 - Lambda - psi_d*c_d*M^(-gamma_d)*(M/m_c)^(1/eta) - lambda_c;

1 - Lambda*(1 + r) - psi_d*theta*c_d*M^(-gamma_d)*(M/m_c)^(1/eta) - lambda_d;

c_d - r*m_d;

y - z_r*a_r^alpha - z_lr*a_lr^alpha;

mu_r + (rho-1)*z_r;

mu_lr + (rho-1)*z_lr;

pi - y - r*m_d;

d - pi;

D - (d/e);

a_r + a_lr - e - m_d;

p*e - xi*(a_r+delta*a_lr);

-p - phi + Lambda*(alpha*z_r*a_r^(alpha - 1) - r)+mu*(p - xi) + lambda_ar;

-p - phi*Lambda*(alpha*z_lr*a_lr^(alpha - 1) - r) - mu*(p - xi*delta) + lambda_ar;

p+phi - Lambda*(1+r) - mu*p;

lambda_d*m_d=0;

lambda_c*m_c=0;

lambda_ar*a_r = 0;

lambda_alr*a_alr = 0;




%% end own model equations

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

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