
%function [ys,check]=BGG_normal_dist_steadystate(ys,exe);
function [ys,params,check]=BGG_normal_dist_steadystate(ys,exo,M_,options_);
%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.
%%


Hbar0=0.5;
K_0=5.1; %4.0;
sigmaE_0=0.98; %0.9;
sigmapsi_0=0.31; %0.3;
psi_0=0.5;
muFF_0=0.015; %0.1;
xiE_0=0.006; %0.1;

x0=[log(Hbar0);  log(K_0); ...
   log(sigmaE_0/(1-sigmaE_0)) ;  log(sigmapsi_0); log(psi_0);log(muFF_0/(1-muFF_0)); log(xiE_0/(1-xiE_0))];
options = optimset('TolFun',1e-5,'TolX',1e-5,'MaxIter', 5000, 'MaxFunEvals', 5000, 'Display','on');

[x,~] = fsolve(@ss_fun_BGG_normal_dist_SW,x0,options);
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)));
g=g_bar;

%
%Then given varaibles above set out the rest of the steady state in recursive form in pure Matlab; note some differences with dynare-matlab 
%
PIE=PIE_bar;
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
A=A_bar;


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


PO=((1-xip*PIEtilde^(zzeta-1))/(1-xip))^(1/(1-zzeta));
MC=(1-1/zzeta)*PO*(1-bettag*(1+g)*xip*PIEtilde^zzeta)/(1-bettag*(1+g)*xip*PIEtilde^(zzeta-1));



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    

PWP=MC;

Deltap=((1-xip)*(((1-xip*PIEtilde^(zzeta-1))/(1-xip))^(1/(1-zzeta)))^(-zzeta))/(1-xip*PIEtilde^zzeta);


Q=1;
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;
Wh=W*Invmarkup;
Q=1;



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

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);
%
%corection
%DG=1/(psi*sigmapsi*sqrt(2*(4*atan(1))))*exp(-(log(psi)+.5*(sigmapsi)^2)^2/(2*(sigmapsi)^2));
%should be
%
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;

%
%Other endogenous Variables
%



U=((C*(1-chi/(1+g)))^(1-sigma_c)*exp((sigma_c-1)*H^(1+sigma_l)/(1+sigma_l)) -1)/(1-sigma_c);


CEequiv=((1.01*(C*(1-chi/(1+g)))^(1-sigma_c)*exp((sigma_c-1)*H^(1+sigma_l)/(1+sigma_l))-1)/(1-sigma_c)-U)/(1-bettag*(1+g));

Omega=U/(1-bettag*(1+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;
%RK=R;
%
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;
if wage>0
JJw=Hd/(1-bettag*xiw*PIEwtilde^muw/PIEwptilde);
Jw=(1/(1-(1/muw)))*Wh*Hd/(1-bettag*xiw*(1+g)*PIEwtilde^muw);
else
JJw=Hd;
Jw=Hd*Wh;   
end

MCS=1;
MRSS=1;
MPS=1;
RPS=1;
IS=1;
Rn=1/bettag*PIE;
ER=1/bettag;

if wage ==0
JJw=Hd;    
Jw=Hd*Wh*MRSS;
W=Wh*MRSS;
Deltaw=1;
WO=Jw/(W*JJw);
Invmarkup=Wh/W;
end


rhoRex=rho*R;

spread=RK-R;
epsAtrendstate=0;

%

%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state vZNues of the endogenous variables of
%% the model.
%%
                                                        % End of the loop.
%%

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end  
