function [ys,params,check] = bank_steadystate(ys,exo,M_,options_)
% function [ys,check] = bank_steadystate(ys,exo)
% computes the steady state for the bank.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)


% 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

% the steady state computation follows 
PI=1.15;
A=1;
or=1;
g=1;
tao=1;
mu1=PI;
i=(PI/betta)-1;
rk=(1/betta)-(1-delta);
rd=PI*(rk+(1-delta))-1;
rl=((1+rd)/alphab*(1-Lambda+roa))+(tao*(1-alphab)/alphab)-1;
mc=(theta-1)/theta;
w=(1-alppha)*(mc*(alppha/(rk*(1+Lambda)))^alppha)^(1/(1-alppha));
c=w/psin;
m=(c*psim)*((1+i)/i);
k=c/((((1-alppha)*(1+Lambda)*rk)/(alppha*w))^(1-alppha)-delta);
n=((1-alppha)*(1+Lambda)*rk*k)/(alppha*w);
x=delta*k;
y=(n^(1-alppha))*(k^alppha);
l=(1+Lambda)*rk*k+w*n;
as=l/(1-Lambda+roa);
PIbar=PI;
mubar=mu1;
ybar=y;
rlbar=rl;
gy_obs=1;
gor_obs=1; 
gg_obs=1;
gc_obs=1;
gk_obs=1;

%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ',1) = ' 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
