%compute the steady state of NK Model 

function [ys,params,check]=NK_SW_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.
%%

Hbar=0.5;

x0=[log(Hbar)];
%x0=[0; 0 ];
options = optimset('TolFun',1e-9,'TolX',1e-9,'MaxIter', 5000, 'MaxFunEvals', 5000, 'Display','on');
[x,~] = fsolve(@ss_fun_NK_SW,x0,options);
H=exp(x(1));


g=g_bar;
PIE=PIE_bar;
A=A_bar;
%H=H_bar;
RPS=RPS_bar;
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


%%%%%%%%%%%%%%%
bettag=betta*(1+g)^(-sigma_c);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R=1.0/(RPS*bettag);
RK=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 Wages
%
WO=1;
Deltaw=1;
Invmarkup=1;   
end   


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

Wh=W*Invmarkup;

C=CY*Y;


%
%
%
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;
Q=1;
%
S=0;
Sdash=0;

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;

IS=1;
Rn=PIE/(RPS*bettag);
ER=1/(RPS*bettag);

spread=RK-R;
N=0.00001*Y;
epsAtrendstate=0;




%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
  %% Here we define the steady state values of the endogenous variables of
%% the model.
%%
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    

