
function F = ss_fun_BGG_normal(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 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Variables to solve in fsolve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H=exp(x(1));
K=exp(x(2));
sigmaE=exp(x(3))/(1+exp(x(3)));
sigmapsi=exp(x(4));
psi=exp(x(5));
muFF=exp(x(6))/(1+exp(x(6))); 
xiE=exp(x(7))/(1+exp(x(7)));

PIE=PIE_bar;
A=A_bar;
g=g_bar;
bettag=betta*(1+g)^(-sigma_c);
R=1.0/bettag;

PIEw=PIE_bar*(1+g);
if index==1
%
%with complete indexing
%
PIEtilde=1;
PIEwtilde=1;
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;
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 Eages
%
WO=1;
Deltaw=1;
Invmarkup=1;   
end 


Hd=H/Deltaw;
YW=A*Hd^(alp)*(K/(1+g))^(1-alp);
KYW=K/YW;

IY=(delta+g)*KYW/(1+g)*Deltap;

Y=YW/Deltap;

I=IY*Y;
G=gy*Y;
W=alp*PWP*YW/Hd;
Q=1;



RK=(1-alp)*PWP*YW*(1+g)/K+1-delta;


p=logncdf(psi,-0.5*(sigmapsi)^2,sigmapsi);
fnG=1- normcdf((0.5*(sigmapsi)^2 -log(psi))/sigmapsi,0,1);
fnGam=fnG+psi*(1-p);
DGam=(1-p);
%
%
DG=1/(sigmapsi*sqrt(2*(4*atan(1))))*exp(-(log(psi)+.5*(sigmapsi)^2)^2/(2*(sigmapsi)^2));


rho=DGam/((fnGam-muFF*fnG)*DGam+(1-fnGam)*(DGam-muFF*DG));
N=(sigmaE+xiE)*(1-fnGam)*RK*Q*K/(1+g);
phi=Q*K/N;
spread=RK-R;
CE=(1-xiE)*(1-sigmaE)*(1-fnGam)*RK*Q*K/(1+g);
RL = PIE*psi*(RK*K/(1+g))/(K/(1+g)-N);
C=Y-CE-G-I-muFF*fnG*RK*Q*K/(1+g);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CY=C/Y;
tax=G;





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


LAMBDA=bettag;
S=0;
Sdash=0;
Q=1;
X=1+g;

JJp=Y/(1-bettag*(1+g)*xip*PIEtilde^(zzeta-1));
JJtilde=(PIEtilde^(zzeta-1))*JJp;
Jp=(1/(1-1/zzeta))*Y*MC/(1-bettag*(1+g)*xip*PIEtilde^(zzeta));
Jtilde=PIEtilde^zzeta*Jp;


MCS=1;
IS=1;
Rn=1/bettag*PIE;
ER=1/bettag;
Wh=W*Invmarkup;


F = [-UH/lam-Wh; rho-rhocalib;  p-pcalib; phi-phicalib; CE/Y-CEcalib;
   RK-rho*R; phi*RK*(fnGam-muFF*fnG)-R*(phi-1);  ]; 

