function [ys,check] = BAINHAP_steadystate(ys,exe)
global M_ lgy_ options_

if isfield(M_,'param_nbr') == 1
NumberOfParameters = M_.param_nbr;
for i = 1:NumberOfParameters
  paramname = deblank(M_.param_names(i,:));
  eval([ paramname ' = M_.params(' int2str(i) ');']);
end
check = 0;
end

beta = 0.99;
sigma = 1;
gamma = 0.5;
phi = 2;
alpha = 0.75;
eta = 0.35;
lamda = 11;
delta = 0.025;
omega = 0.01;
alpha_F = 0.1;
psi_I = 12;
psi_D = 0.0075;
psi_H = 120;
psi_X = 120;
psi_M = 120;
epsilon_pi = 1.5;
rho = 0.5;
chi_D = 0.9;
mu_D = 0.2;
theta_D = 0.9873;
chi_F = 0.5;
mu_F = 0.12;
theta_F = 0.9966;

xi = 1 ;


p_D_ss = 1;
p_F_ss = 1; 
 
C_D_ss = 1;
C_F_ss = 1;

I_D_ss = 1; 
I_F_ss =1;

i_D_ss = 1/beta -1; 
i_F_ss = 1/beta -1;

pi_D_ss = 1;
pi_F_ss = 1;

K_D_ss = I_D_ss/delta;
K_F_ss = I_F_ss/delta;

Q_D_ss = p_D_ss;
Q_F_ss = p_F_ss;

PHI_D_ss = 0.02;
PHI_F_ss = 0.005;

r_KD_ss = (1+i_F_ss)*(1+PHI_D_ss);
r_KF_ss = (1+i_D_ss)*(1+PHI_F_ss);

r_D_ss = r_KD_ss-(1-delta);
r_F_ss = r_KF_ss-(1-delta);
 

THETA_D_ss = 1;
THETA_F_ss = 1;

h_D_ss = exp(log((1-eta)*(1-omega)*r_D_ss*K_D_ss/eta)/(phi+1));
h_F_ss = exp(log((1-eta)*(1-omega)*r_F_ss*K_F_ss/eta)/(phi+1));
 
 
w_D_ss = h_D_ss^phi;
w_F_ss = h_F_ss^phi;

n_D_ss = h_D_ss^(1-omega);
n_F_ss = h_F_ss^(1-omega);
 
y_D_ss = n_D_ss^(1-eta)*K_D_ss^eta;
y_F_ss = n_F_ss^(1-eta)*K_F_ss^eta;
 
mc_D_ss = r_D_ss*K_D_ss/(eta*y_D_ss);
mc_F_ss = r_F_ss*K_F_ss/(eta*y_F_ss);

w_ED_ss = ((1 - eta)*omega*y_D_ss * mc_D_ss);
w_EF_ss = ((1 - eta)*omega*y_F_ss * mc_F_ss);

p_HD_ss = lamda/(lamda - 1)*mc_D_ss;
p_HF_ss = lamda/(lamda - 1)*mc_F_ss;

C_HD_ss = (alpha)*(p_HD_ss/p_D_ss)^(-gamma);
C_HF_ss =(alpha_F)*(p_HF_ss/p_F_ss)^(-gamma) ;

C_MD_ss = exp(log((1-(alpha)^(1/gamma)*C_HD_ss^((gamma-1)/gamma))/((1-alpha)^(1/gamma)))*(gamma/(gamma-1)));
C_MF_ss = exp(log((1-(alpha_F)^(1/gamma)*C_HF_ss^((gamma-1)/gamma))/((1-alpha_F)^(1/gamma)))*(gamma/(gamma-1)));


p_MD_ss = exp(log(C_MD_ss/(1-alpha))/(-gamma));
p_MF_ss = exp(log(C_MF_ss/(1-alpha))/(-gamma));

I_HF_ss = (alpha_F)*(p_HF_ss/p_F_ss);

I_MF_ss = (1-alpha_F)*(p_MF_ss/p_F_ss);

p_XD_ss = lamda/(lamda - 1)*mc_D_ss;
p_XF_ss = lamda/(lamda - 1)*mc_F_ss;

S = 0.1207;
RHO_D_ss = 1;
om_D_ss = 0.603892;
z_D_ss = 1 - normcdf(((log(om_D_ss)+0.5*S^2)/S) - S)- om_D_ss*(1-normcdf(((log(om_D_ss)+0.5*S^2)/S))) ;
g_D_ss = RHO_D_ss*om_D_ss/RHO_D_ss*(1-normcdf(((log(om_D_ss/RHO_D_ss)+0.5*S^2)/S))) + (1-mu_D)*normcdf(((log(om_D_ss/RHO_D_ss)+0.5*S^2)/S) - S);
zph_D_ss = -(1 - normcdf(((log(om_D_ss)+0.5*S^2)/S)));
gph_D_ss = RHO_D_ss*(1-normcdf(((log(om_D_ss/RHO_D_ss)+0.5*S^2)/S)) - mu_D/(sqrt(2*3.141592654)*S)*exp(-((log(om_D_ss/RHO_D_ss)+0.5*S^2)/S)^2/2)); 

S_F = 0.12722440917;
RHO_F_ss = 1;
om_F_ss = 0.640610102;
z_F_ss = 1 - normcdf(((log(om_F_ss)+0.5*S_F^2)/S_F) - S_F)- om_F_ss*(1-normcdf(((log(om_F_ss)+0.5*S_F^2)/S_F))) ;
g_F_ss = RHO_F_ss*om_F_ss/RHO_F_ss*(1-normcdf(((log(om_F_ss/RHO_F_ss)+0.5*S_F^2)/S_F))) + (1-mu_F)*normcdf(((log(om_F_ss/RHO_F_ss)+0.5*S_F^2)/S_F) - S_F);
zph_F_ss = -(1 - normcdf(((log(om_F_ss)+0.5*S_F^2)/S_F)));
gph_F_ss = RHO_F_ss*(1-normcdf(((log(om_F_ss/RHO_F_ss)+0.5*S_F^2)/S_F)) - mu_F/(sqrt(2*3.141592654)*S_F)*exp(-((log(om_F_ss/RHO_F_ss)+0.5*S_F^2)/S_F)^2/2)); 


D_ED_ss = (r_KD_ss*K_D_ss*g_D_ss)/(1+i_F_ss) ;
D_EF_ss = (r_KF_ss*K_F_ss*g_F_ss)/(1+i_D_ss) ;

nw_D_ss = K_D_ss - D_ED_ss;
nw_F_ss = K_F_ss -D_EF_ss;

C_ED_ss = 1;
C_EF_ss = 1;

C_EHD_ss = C_HD_ss*C_ED_ss;
C_EHF_ss = C_HF_ss*C_EF_ss;

C_EMF_ss = C_MF_ss*C_EF_ss;

y_XD_ss = alpha_F*(p_XD_ss)^(-gamma)*y_F_ss;
y_XF_ss = alpha*(p_XF_ss)^(-gamma)*y_D_ss;

y_HD_ss = y_D_ss - y_XD_ss;
y_HF_ss = y_F_ss - y_XF_ss;

y_MD_ss = C_MD_ss + C_MD_ss*C_ED_ss + ((1-alpha)*(p_MD_ss/p_D_ss)) + (1-alpha)*p_MD_ss^(-gamma)*(1-z_D_ss-g_D_ss)*r_KD_ss*Q_D_ss*K_D_ss;
y_MF_ss = C_MF_ss + C_EMF_ss + I_MF_ss + (1-alpha_F)*p_MF_ss^(-gamma)*(1-z_F_ss-g_F_ss)*r_KF_ss*Q_F_ss*K_F_ss;
 
D_HD_ss =((p_XD_ss*y_XD_ss - p_MF_ss*y_MD_ss)- i_F_ss*D_ED_ss)/((1+i_F_ss)*psi_D -1);
D_HF_ss =((p_XF_ss*y_XF_ss - p_MD_ss*y_MF_ss)- i_D_ss*D_EF_ss)/((1+i_D_ss)*psi_D -1);
 

                                                          
% save steady state values of variables


C_HD    =   C_HD_ss ;
C_MD    =   C_MD_ss ;
C_D     =   C_D_ss ;
C_HF    =   C_HF_ss ;
C_MF    =   C_MF_ss ;
C_F     =   C_F_ss ;
p_HD    =   p_HD_ss ;
p_MD    =   p_MD_ss ;
p_D     =   p_D_ss ;
        
p_HF    =   p_HF_ss ;
p_MF    =   p_MF_ss ;
p_F     =   p_F_ss ;
        
        
        
mc_D    =   mc_D_ss ;
mc_F    =   mc_F_ss ;
        

p_XD    =   p_XD_ss ;
p_XF    =   p_XF_ss ;
        
I_D     =   I_D_ss ;
        
I_F     =   I_F_ss ;
I_HF    =   I_HF_ss ;
I_MF    =   I_MF_ss ;
        
K_D     =   K_D_ss ;
K_F     =   K_F_ss ;
        
Q_D     =   Q_D_ss ;
Q_F     =   Q_F_ss ;
        
THETA_D     =   THETA_D_ss ;
THETA_F     =   THETA_F_ss ;
        
RHO_D   =   RHO_D_ss ;
RHO_F   =   RHO_F_ss;
om_D    =   om_D_ss ;
om_F    =   om_F_ss ;
z_D     =   z_D_ss   ;
g_D =   g_D_ss;
zph_D   =   zph_D_ss  ;
gph_D   =   gph_D_ss ; 
z_F     =   z_F_ss  ;
g_F     =   g_F_ss ;
zph_F   =   zph_F_ss ;
gph_F   =   gph_F_ss ; 
        
PHI_D   =   PHI_D_ss ;
PHI_F   =   PHI_F_ss ;
        
i_D     =   i_D_ss ; 
i_F     =   i_F_ss ;
        
r_KD    =   r_KD_ss ;
r_KF    =   r_KF_ss ;
        
D_ED    =   D_ED_ss ;
D_EF    =   D_EF_ss;
        
C_ED    =   C_ED_ss ;
C_EF    =   C_EF_ss ;
        
C_EHD   =   C_EHD_ss ;
C_EHF   =   C_EHF_ss ;
C_EMF   =   C_EMF_ss ;
        
y_D     =   y_D_ss ;
y_HD    =   y_HD_ss ;
y_XD    =   y_XD_ss ;
        
y_F     =   y_F_ss ;
y_HF    =   y_HF_ss ;
y_XF    =   y_XF_ss  ;
        
y_MD    =   y_MD_ss ;
y_MF    =   y_MF_ss ;
        
w_ED    =   w_ED_ss ;
w_EF    =   w_EF_ss ;
        
nw_D    =   nw_D_ss;
nw_F    =   nw_F_ss;
        
       
w_D     =   w_D_ss  ;
w_F     =   w_F_ss  ;
        
n_D     =   n_D_ss ;
n_F     =   n_F_ss ;
        
h_D     =   h_D_ss ;
h_F     =   h_F_ss ;
        
D_HD    =   D_HD_ss  ;
D_HF    =   D_HF_ss ;
        
        
pi_D    =   pi_D_ss  ;
pi_F    =   pi_F_ss  ;
        
r_D     =   r_D_ss ;
r_F     =   r_F_ss ;
        
psi_Dt = psi_D;
 

%%
for iter = 1:length(M_.params)
  eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

if isfield(M_,'param_nbr') == 1

if isfield(M_,'orig_endo_nbr') == 1
NumberOfEndogenousVariables = M_.orig_endo_nbr;
else
NumberOfEndogenousVariables = M_.endo_nbr;
end
ys = zeros(NumberOfEndogenousVariables,1);
if options_.linear
    
else
for i = 1:NumberOfEndogenousVariables
  varname = deblank(M_.endo_names(i,:));
  eval(['ys(' int2str(i) ') = ' varname ';']);
end
end

else
ys=zeros(length(lgy_),1);
for i = 1:length(lgy_)
    ys(i) = eval(lgy_(i,:));
end
check = 0;
end