function [ys,check] = extmodelsp_steadystate(ys,exe)

%% Load model parameters
global M_
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end         
check = 0;
%%

R_h = 1/beta_h;

lambda_b = (1-beta_b*R_h)/(1-beta_b*rho_d);


R_e = (1/beta_b)*(1-((1-beta_b)*rho_d+(1-rho_d)*gamma_e)*lambda_b);

R_s = (1/beta_b)*(1-((1-beta_b)*rho_d+(1-rho_d)*gamma_s)*lambda_b);

lambda_e = (1-beta_e*R_e)/(1-beta_e*rho_e);

lambda_s = (1-beta_s*R_s)/(1-beta_s*rho_s);

R_k = 1/beta_e*(1-lambda_e*(1-rho_e)*m_k) - (1-delta_ke_ss);   

R_v = 1/beta_e*(1-lambda_e*m_h/R_e*(1-rho_e)) - 1;

R_m = 1/beta_h - (1-delta_kh_ss);   



oo1 = j/(1-beta_h);
oo2 = j/(1-beta_s-lambda_s*(1-rho_s)*m_s/R_s);
oo3 = 1/(1+(1-1/R_s)*m_s*oo2);
oo4 = gamma_e*(R_h-1)*(m_h*nu/(R_e*R_v)+m_k*alp*mu/(R_k)-m_n*(1-alp-nu)/(1+m_n*lambda_e));     
oo5 = gamma_s*m_s/R_s*(oo2*(R_h-1))*oo3;
oo6 = R_m - delta_ke_ss;   
oo7 = (1-alp-nu)*(1-sigma)/((1+m_n*lambda_e));
oo8 = (1-mu)*alp/(R_m);


N_s = 1/(1+oo3*tau);   

zeta1 = oo6*oo8/oo7 + 1 + oo4/oo7 + oo5*sigma/(1-sigma);        

N_h = 1/(1+tau*zeta1); 

ctoys = oo3*(1-alp-nu)/(1+m_n*lambda_e)*sigma;
ctoyh = zeta1*oo7 ;
iktoye = delta_ke_ss*alp*mu/(R_k);    
iktoyh = delta_ke_ss*alp*(1-mu)/(R_m);
denh = oo1*ctoyh + oo2*ctoys + nu/(R_v);

H_h = oo1*ctoyh / denh ;
H_s = oo2*ctoys / denh ;
H_e = (nu/(R_v)) / denh ;

Y = (alp*(1-mu)/(R_m))^(alp*(1-mu)/(1-alp))*...
    (alp*mu/(R_k))^(alp*mu/(1-alp))*...
    H_e^(nu/(1-alp))*...
    (N_h^(1-sigma)*N_s^sigma)^((1-alp-nu)/(1-alp));

K_h = alp*(1-mu)*Y/(R_m);
K_e = alp*mu*Y/(R_k);

q = nu*Y/(R_v)/H_e;

L_e = m_h*q/R_e*H_e + m_k*K_e - m_n*(1-alp-nu)/(1+m_n*lambda_e)*Y;  

W_h = ((1-sigma)*(1-alp-nu)/(1+m_n*lambda_e))*Y/N_h;
W_s = (sigma*(1-alp-nu)/(1+m_n*lambda_e))*Y/N_s;

L_s = m_s/R_s*oo2*oo3*W_s*N_s;

D = gamma_e*L_e + gamma_s*L_s;

C_h = (R_m-delta_ke_ss)*K_h + (R_h-1)*D + W_h*N_h;       

C_b = (R_e-1)*L_e + (R_s-1)*L_s - (R_h-1)*D;

C_e = Y - delta_ke_ss*K_e - (R_e-1)*L_e - W_h*N_h - W_s*N_s - R_m*K_h;

C_s = oo3*W_s*N_s;

z_kh = 1;

z_ke = 1;

eps_e = 0;
eps_h = 0;

%A_j = 1;
%A_k = 1;
%A_me = 1;
%A_mh = 1;
%A_p = 1;
%A_z = 1;




%% DO NOT CHANGE THIS PART.
%%
%% Store model parameters
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval(['M_.params(' int2str(i) ') = ' paramname ';']);       %    Get the value of parameter i.
end   


%% Here we define the steady state values of the endogenous variables of
%% the model.
%%
NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for i = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
  eval(['ys(' int2str(i) ') = ' varname ';']);                %    Get the steady state vZNue of this variable.
end                                                           % End of the loop.
%%
