function [ys,check] = efdc_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

 

betaHUU=0.995;
betaEUU=0.98;
betaGUU=0.98;
aHUU=0.63;
aEUU=0.63;
jUU=0.06;
deltaUU=0.04;
mukUU=0.5;
mulUU=0.1;
mugUU=0.02; 
mEUU=0.3;
mGUU=0.5;
epsUU=11;
nuUU=0.2;
omega2UU=0.07;
etaUU=0.2;
 


piUU=1;
AUU=1;
qkUU=1;
epscUU=1;
epshUU=1;
epslUU=1;
xUU=epsUU/(epsUU-1);
RUU=1/betaHUU;
zetaH=(1-betaHUU*aHUU)/(1-aHUU);
zetaE=(1-betaEUU*aEUU)/(1-aEUU);
zeta1=(1/betaEUU-1+deltaUU-((1-betaEUU*RUU)*mEUU*(1-deltaUU))/(betaEUU*RUU))^(-1)*mukUU/xUU;
zeta2=betaEUU*mulUU/xUU*(1-betaEUU-mEUU*(1-betaEUU*RUU)/RUU)^(-1);
zeta3=jUU/(zetaH*(1-betaHUU));
zeta4=mEUU/RUU*((1-deltaUU)*zeta1+zeta2);
zeta5=betaGUU*mugUU/(xUU*(1-betaGUU*(1-deltaUU)));
zeta6=((1-betaGUU)*RUU-(1-betaGUU*RUU)*mGUU)/(omega2UU*RUU);
zeta7=(mugUU/xUU+etaUU-deltaUU*zeta5)/((1-omega2UU)*zeta6+(RUU-1)*mGUU/RUU);
zeta8=mGUU/RUU*zeta7;
zeta9=(1-mukUU-mulUU-mugUU)/xUU;
zeta10=zeta9+(RUU-1)*(zeta4+zeta8)+(1-1/xUU)-etaUU;
zeta11=zeta3*zeta10+zeta2+zeta7;
zeta12=zeta9*zetaH/zeta10;
zeta13=nuUU*zeta6*zeta7/deltaUU+zeta5;
nUU=zeta12^(1/(1+phiUU));
zeta14=zeta1^mukUU*(zeta2/zeta11)^mulUU*zeta13^mugUU*nUU^(1-mukUU-mulUU-mugUU);
zeta15=1/xUU-zeta9-(RUU-1)*zeta4-deltaUU*zeta1-mugUU/xUU;
lrUU=zeta7/zeta11;
lsUU=1-lrUU;
leUU=zeta2/zeta11;
yUU=zeta14^(1/(1-mukUU-mugUU));
qlUU=zeta11*yUU;
lhUU=lsUU-leUU;
keUU=zeta1*yUU;
iUU=deltaUU*keUU;
bEUU=zeta4*yUU;
bGUU=zeta8*yUU;
bUU=bEUU+bGUU;
dUU=bUU;
cHUU=zeta10*yUU;
wUU=zeta9*yUU/nUU;
kgUU=zeta13*yUU;
cEUU=zeta15*yUU;
cUU=cHUU+cEUU;
taxUU=etaUU*yUU;
lamHUU=zetaH/cHUU;
lamEUU=zetaE/cEUU;
thEUU=(1-betaEUU*RUU)*lamEUU/RUU;
FUU= (1-1/xUU)*yUU;
lamGUU=1/(zeta6*qlUU*lrUU);
gcUU=(1-nuUU-omega2UU)/lamGUU;
giUU=deltaUU*kgUU;
chiGUU=lamGUU-nuUU/giUU;
rgUU= mugUU*yUU/(xUU*kgUU);
thGUU=lamGUU*(1-betaGUU*RUU)/RUU;
PHIUU=lamHUU*yUU/(1-thetaUU*betaHUU);
OMEGAUU=PHIUU;


lamHU =log(lamHUU);
cHU =log(cHUU);
wU =log(wUU);
nU =log(nUU);
RU =log(RUU);
dU =log(dUU);
FU =log(FUU);                                 
lamEU =log(lamEUU);
lamGU =log(lamGUU);
cEU =log(cEUU);
thEU =log(thEUU);
yU =log(yUU);
xU =log(xUU);
keU =log(keUU); 
bEU =log(bEUU);
bGU = log(bGUU);
thGU=log(thGUU);
iU =log(iUU);                                                                                                       
gcU =log(gcUU);   
kgU=log(kgUU);
giU=log(giUU);
rgU=log(rgUU);
chiGU=log(chiGUU);
cU =log(cUU);
bU =log(bUU);  
piU =log(piUU);
leU =log(leUU);
lhU =log(lhUU);
lsU =log(lsUU);
lrU =log(lrUU);
AU =log(AUU);
qkU =log(qkUU);
qlU=log(qlUU);
epscU=log(epscUU);
epshU=log(epshUU);
epslU=log(epslUU);
PHIU=log(PHIUU);
OMEGAU=log(OMEGAUU);
taxU=log(taxUU);
piU_ob=log(piUU);
cU_ob=0;
iU_ob=0;
yU_ob=0;
RU_ob=0;
qlU_ob=0;

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

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