function[ys,check] = FO_OE_ss_steadystate(ys,exo)

% computes the steady state for the FO_OE_ss.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 for the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose 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 

% the steady state computation

%%% here I use the first solution for dt - with plus

dt=((((1-alfa*csi)/(alfa*mu))^2-(1-alfa*csi)/(alfa*mu)+1)^(1/2)...
    +2*(1-alfa*csi)/(alfa*mu)-1)/(6*(phi+gama)/mu);

%%%%% second solution for dt is
%%%% dt=(-(((1-alfa*csi)/(alfa*mu))^2-(1-alfa*csi)/(alfa*mu)+1)^(1/2)+2*(1-alfa*csi)/(alfa*mu)-1)/(6*(phi+gama)/mu);
 
%%% now, derived equations for d and dstar, as functions of dt
d=-dt*(1/alfa-(phi+gama)*dt-csi)/(1/alfa-2*(phi+gama)*dt-csi);
dstar=dt*(1/alfa-(phi+gama)*dt-csi-mu)/(1/alfa-2*(phi+gama)*dt-csi-mu);

rd = gama*dt;                                              %% functional form returns on deposits
ri = 1/alfa-((phi/2)*dt);                          %% implicit function returns on firms' projects
rl = 1/alfa - phi*dt;                                        %% functional form returns on loans
pr =1- alfa*ri;                               %% functional form probability

%set the limitation condition for pr 
 if pr>=1 || pr <0  % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form 
 check=1; %set failure indicator 
  return; %return without updating steady states 
 end
 
 

pai = pr*((1+rl)-(1+rd)-csi)*d;                     %%% current banks' profits  domestic destination market
paistar = pr*((1+rl)-(1+rd)-csi-mu)*dstar;          %%% current banks' profits  foreign destination market

ne=0.00004;                                   %%%% SET initial value ne
n = ne*(1-pho)/pho;
nstar =(dt*(1-pho)- n*d)/dstar;              
nestar = nstar*pho/(1-pho);
nt = n + nstar;                     %% aggregate banks
net = ne + nestar;               %% aggregate new banks
a=0;
m=0;


%%%%% end model %%%%%
%%% update parameters 
 
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