
function F = ss_fun(x)

    %% Load model parameters
global M_
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
 
end 

H=exp(x(1));
g=g_bar;
PIE=PIE_bar;
A=A_bar;

bettag=betta*(1+g)^(-sigma_c);
RK=1.0/bettag;

PIEw=PIE_bar*(1+g);
if index==1
%
%with complete indexing
%
PIEtilde=1;
PIEwtilde=1+g;
PIEwptilde=1;

else
PIEtilde=PIE^(1-gammap);
PIEwtilde=PIEw/PIE^gammaw;
PIEwptilde=PIE^(1-gammaw);
end

MC=(1-1/zzeta)*(1-xip*bettag*(1+g)*PIEtilde^zzeta)/(1-xip*bettag*(1+g)*PIEtilde^(zzeta-1))*(((1-xip*PIEtilde^(zzeta-1))...
/(1-xip))^(1/(1-zzeta)));
PWP=MC;
KYW=PWP*(1+g)*(1-alp)/(RK-1+delta);
Deltap=((1-xip)*(((1-xip*PIEtilde^(zzeta-1))/(1-xip))^(1/(1-zzeta)))^(-zzeta))/(1-xip*PIEtilde^zzeta);
if wage>0
WO=((1-xiw*PIEwtilde^(muw-1))/(1-xiw))^(1/(1-muw));
Deltaw=(1-xiw)/(1-xiw*PIEwtilde^muw)*WO^(-muw);
Invmarkup=(1-1/muw)*WO*(1-bettag*xiw*(1+g)*PIEwtilde^muw)/(1-bettag*xiw*PIEwtilde^muw*PIEwptilde^(-1));
else
%Sticky Wages
%
WO=1;
Deltaw=1;
Invmarkup=1;   
end


IY=(1-alp)*(delta+g)*(MC*Deltap)/(RK-1+delta);
CY=1-IY-gy;
Hd=H/Deltaw;
YW=A*Hd*(KYW/(1+g))^((1-alp)/alp);

Y=YW/Deltap;

K=KYW*YW;
I=IY*Y;
G=gy*Y;
W=alp*PWP*YW/Hd;
C=CY*Y;

UC=(C*(1-chi)/(1+g))^(-sigma_c)*exp((sigma_c-1)*H^(1+sigma_l)/(1+sigma_l));
UH=-UC*H^(sigma_l)*C*(1-chi/(1+g));
%
if  habit==1
%
%external habit
%
lam=UC;
%
else
%
%internal habit
%
lam=UC*(1-bettag*chi);
end



Wh=W*Invmarkup;
F = [-UH/lam-Wh ] ;

