function [ys,check] = chay1_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 = 2;
phi = 2; 
alpha = 0.315; 
eta = 0.35; 
lamda = 11; 
delta = 0.025;
omega = 0.01; 
alpha_F = 0.035; 
psi_I = 12; 
psi_D = 0; 
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	=	40;
K_F_ss	=	40;
Q_D_ss	=	1;
Q_F_ss 	=	1;
r_KD_ss	=	(1+i_D_ss)*(1+PHI_D_ss);
r_KF_ss 	=(1+i_F_ss)*(1+PHI_F_ss);
r_D_ss 	=	r_KD_ss-(1-delta);
r_F_ss 	=	r_KF_ss-(1-delta);
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));

n_D_ss = h_D_ss^(1-omega);
n_F_ss = h_F_ss^(1-omega);

w_D_ss = h_D_ss^phi;
w_F_ss = h_F_ss^phi;

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_F))/(-gamma));

I_HD_ss 	=	(alpha)*(p_HD_ss/p_D_ss)^(-gamma)*I_D_ss;
I_HF_ss 	=	(alpha_F)*(p_HF_ss/p_F_ss)^(-gamma)*I_F_ss;
I_MD_ss 	=	(1-alpha)*(p_MD_ss/p_D_ss)^(-gamma)*I_D_ss;
I_MF_ss 	=	(1-alpha_F)*(p_MF_ss/p_F_ss)^(-gamma)*I_F_ss;

s_D_ss 	=	p_MD_ss*(lamda-1)/lamda/p_F_ss;
s_F_ss = 1/s_D_ss;	

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


THETA_D_ss = beta;		
THETA_F_ss =beta;


PHI_D_ss =  0.02+7.9431e-06;
PHI_F_ss = 0.005-4.3497e-05;
om_D_ss	=	 0.2837 ;
om_F_ss	=	  0.2837 ;
S	=	0.5872;
S_F	=	0.4885;
RHO_D_ss= 1 ;
RHO_F_ss = 1;

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)); 

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*Q_D_ss*K_D_ss*g_D_ss)/((1+i_F_ss)*s_D_ss) ;
D_EF_ss = (r_KF_ss*Q_F_ss*K_F_ss*g_F_ss)/((1+i_D_ss)*s_F_ss);

nw_D_ss = (theta_D*(r_KD_ss*Q_D_ss*K_D_ss*(z_D_ss+g_D_ss)-(1+i_F_ss)*s_D_ss*D_ED_ss)+w_ED_ss)/p_D_ss;
nw_F_ss = (theta_F*(r_KF_ss*Q_F_ss*K_F_ss*(z_F_ss+g_F_ss)-(1+i_D_ss)*s_F_ss*D_EF_ss)+w_EF_ss)/p_F_ss;

C_ED_ss = (1-theta_D)*(r_KD_ss*Q_D_ss*K_D_ss*(z_D_ss+g_D_ss) - (1+i_F_ss)*s_D_ss*D_ED_ss);
C_EF_ss = (1-theta_F)*(r_KF_ss*Q_F_ss*K_F_ss*(z_F_ss+g_F_ss) - (1+i_D_ss)*s_F_ss*D_EF_ss);

C_EHD_ss	=	C_ED_ss*C_HD_ss;
C_EHF_ss = (alpha_F)*(p_HF_ss/p_F_ss)^(-gamma)*C_EF_ss;
C_EMD_ss	=	C_ED_ss*C_MD_ss;
C_EMF_ss	=(1-alpha_F)*(p_MF_ss/p_F_ss)^(-gamma)*C_EF_ss;


y_HD_ss = C_HD_ss + C_EHD_ss + I_HD_ss + (alpha)*(p_HD_ss/p_D_ss)^(-gamma)*((1-z_D_ss-g_D_ss)*r_KD_ss/p_D_ss*Q_D_ss*K_D_ss);
y_HF_ss = C_HF_ss + C_EHF_ss + I_HF_ss + (alpha_F)*(p_HF_ss/p_F_ss)^(-gamma)*((1-z_F_ss-g_F_ss)*r_KF_ss/p_F_ss*Q_F_ss*K_F_ss);
y_XD_ss	=	y_D_ss - y_HD_ss;
y_XF_ss	=	y_F_ss - y_HF_ss;
y_MD_ss =C_MD_ss + C_EMD_ss + I_MD_ss + (1-alpha)*(p_MD_ss/p_D_ss)^(-gamma)*((1-z_D_ss-g_D_ss)*r_KD_ss/p_D_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/p_F_ss)^(-gamma)*((1-z_F_ss-g_F_ss)*r_KF_ss/p_F_ss*Q_F_ss*K_F_ss);
D_HD_ss 	=	(p_XD_ss*y_XD_ss-p_F_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_D_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_EMD   =   C_EMD_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;
I_HD = I_HD_ss;
I_MD = I_MD_ss;
s_D = s_D_ss;
s_F = s_F_ss;

%%
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