function [ys,params,check] = model_nosb_test_steadystate(ys,exo,M_,options_)

% 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

aux1p = (1-betap*ap)/(1-ap); %SS relationship between cp and lambdap
aux1e = (1-betae*ae)/(1-ae); %SS relationship between ce and lambdae

pie = 1;
Lambdap = betap;                                            
Rcd = 1/betap-1;                                              
R = Rcd + 0.0035/4; %35 annual basis points spread (GM 2023)
epsiloncd = Rcd/(Rcd-R);                                       
z=1;
mui = 1;

% commercial banking
Rcb = R;                                                       
Rce = R + 0.0240/4; %240 annual basis points spread (GM 2023)
epsilonce = Rce/(Rce-R);                              

% labour market 
l = 1; 
vw = 1;                                                       
ld = l; 

%retailers
epsilonp = epsilonp_ss;
pstar = 1;                                                   
pe = (epsilonp-1)/epsilonp;                                 
vp = 1;

%capital producers
q = 1;                                                          

%entrepreneurs
A = 1; 

%derivation of capital to gdp ratio 
mc = mc_ss;
kovery = ((((1-(1-deltak)*(mc*(1-betae*(1+Rce))/(1+Rce)))/betae)-(1-deltak))/(alphae*pe))^(-1);
Y = ((kovery^alphae)*(ld^(1-alphae)))^(1/(1-alphae)); 
ye = Y;
k = Y*kovery; %from capital to gdp ratio
w = (1-alphae)*pe*Y/l; 

wstar = w;                                                      
epsilonw = epsilonw_ss;
mrs = ((epsilonw-1)/epsilonw)*w;                              
f1 = (mrs*ld*(w^epsilonw))/(1-phiw*Lambdap);                 
f2 = (ld*(w^epsilonw))/(1-phiw*Lambdap);                     
x1 = (pe*Y)/(1-phip*Lambdap);                                
x2 = (Y)/(1-phip*Lambdap); 

bce = (mc*(1-deltak)*k)/(1+Rce);                                  

ce = pe*ye - deltak*k - w*ld - Rce*bce;       
lambdae = aux1e*(ce^(-1));                                   
deltace = ((1-betae*(1+Rce))/(1+Rce))*lambdae;                     

I = deltak*k;                                                    

nuc = nuc_ss;
kc = nuc*bce; %from capital to asset ratio                                       
dc = bce - kc;                                       
jc = Rce*bce - Rcd*dc;                  
sigmac = jc/kc;                                             

C = Y - I - jc;                                              
cp = C-ce;                                                    
lambdap = aux1p*(cp^(-1));                                 
chi = lambdap*mrs/(l^etap); 

% welfare
up = log(cp-ap*cp) - chi*(l^(1+etap))/(1+etap);
ue = log(ce-ae*ce);
valp = up/(1-betap);
vale = ue/(1-betae);
val = (1-betae)*vale + (1-betap)*valp;


%%
%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state values of the endogenous variables of
%% the model.
%%
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
