function [ys,check] = hm_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;

piUU=1;
epscUU=1;
epshUU=1;
epslUU=1;
betaPUU=0.995;
betaIUU=0.98;
betaEUU=0.98;
betaHUU=0.98;
betaGUU=0.98;
aPUU=0.63;
aIUU=0.63;
aEUU=0.63;
aHUU=0.63;
jUU=0.2;
zetaUU=1.5;   
etaUU=0.5;
omegaUU=0.75;  
deltaUU=0.04;
deltasUU=0.025;
deltahUU=0.025;
muckUU=0.55;
muclUU=0.05;
mugUU=0.07; 
muhkUU=0.4;
gammaUU=0.85;
mIUU=0.7;
mEUU=0.1;
mHUU=0.5;
epsUU=11; 
tauUU=0.2;
nuUU=0.4;
tauhUU=0;


%{

%}

AcUU=1;
AhUU=1;
qckUU=1;
qhkUU=1;
epscUU=1;
epshUU=1;
xUU=epsUU/(epsUU-1);
RUU=1/betaPUU;
zetaP=(1-betaPUU*aPUU)/(1-aPUU);
zetaI=(1-betaIUU*aIUU)/(1-aIUU);
zetaE=(1-betaEUU*aEUU)/(1-aEUU);
zetaH=(1-betaHUU*aHUU)/(1-aHUU);
zeta1=(1/betaEUU-1+deltaUU-((1-betaEUU*RUU)*mEUU*(1-deltaUU))/(betaEUU*RUU))^(-1)*muckUU/xUU;
zeta2=betaEUU*muclUU/xUU*(1-betaEUU-mEUU*(1-betaEUU*RUU)/RUU)^(-1);
zeta3=muhkUU*deltasUU*gammaUU*(1-betaHUU*(1-deltahUU))/(1-betaHUU*(1-deltasUU))*(1/betaHUU-1+deltaUU-((1-betaHUU*RUU)*mHUU*(1-deltaUU))/(betaHUU*RUU))^(-1);
zeta4=jUU/(zetaP*(1+tauhUU-betaPUU*(1-deltahUU)));
zetaR=(zetaUU*(1+tauhUU-betaIUU*(1-deltahUU)-mIUU*(1-betaIUU*RUU)/RUU)/(1+tauhUU-betaPUU*(1-deltahUU)))^(1/etaUU)*(1-omegaUU)/(zetaUU*omegaUU);
zeta5=jUU*omegaUU^etaUU/(zetaI*(omegaUU^etaUU+(1-omegaUU)^etaUU*(zetaUU*zetaR)^(1-etaUU))*(1+tauhUU-betaIUU*(1-deltahUU)-mIUU*(1-betaIUU*RUU)/RUU));
zeta6=mEUU/RUU*((1-deltaUU)*zeta1+zeta2);
zeta7=mHUU/RUU*(1-deltaUU)*zeta3;
zeta8=(1-muckUU-muclUU-mugUU)*alphaUU/xUU;
zeta9=(1-muckUU-muclUU-mugUU)*(1-alphaUU)/xUU;
zeta10=alphaUU*(1-muhkUU)*deltasUU*gammaUU*(1-betaHUU*(1-deltahUU))/(1-betaHUU*(1-deltasUU));
zeta11=(1-alphaUU)*(1-muhkUU)*deltasUU*gammaUU*(1-betaHUU*(1-deltahUU))/(1-betaHUU*(1-deltasUU));
zeta12=(1+tauhUU-betaPUU*(1-deltahUU))/zetaUU;
zeta13=1+(tauhUU+deltahUU)*zeta4-(zeta10-tauUU*deltahUU+zeta7*(RUU-1))*zeta4;
zeta14=(tauhUU+deltahUU)*zeta5*zetaR-zetaUU*zeta12*zetaR*zeta5-mIUU*(RUU-1)*zeta5/RUU-(zeta10-tauUU*deltahUU+zeta7*(RUU-1))*(1+zetaR)*zeta5;
zeta15=zeta8+1-1/xUU+zeta6*(RUU-1)-tauUU;
zeta16=-zeta4*zeta11;
zeta17=1+(tauhUU+deltahUU+(RUU-1)*mIUU/RUU)*zeta5+zetaUU*zeta12*zetaR*zeta5-(1+zetaR)*zeta5*zeta11;
zeta18=(zeta15*zeta17-zeta9*zeta14)/(zeta13*zeta17-zeta14*zeta16);
zeta19=(zeta9*zeta13-zeta15*zeta16)/(zeta13*zeta17-zeta14*zeta16);
zeta20=zeta4*zeta18+(1+zetaR)*zeta5*zeta19;
zeta21=(zeta10*zeta20/zeta8)^(1/(1+xiPUU));
zeta22=(zeta8*zetaP/zeta18*(1+zeta21^(1+xiPUU))^((xiPUU-phiPUU)/(1+xiPUU)))^(1/(1+phiPUU));
zeta23=(zeta11*zeta20/zeta9)^(1/(1+xiIUU));
zeta24=(zeta9*zetaI/zeta19*(1+zeta23^(1+xiIUU))^((xiIUU-phiIUU)/(1+xiIUU)))^(1/(1+phiIUU));
zeta25=(1-gammaUU)*(1-betaHUU*(1-deltahUU))/(1-betaHUU);
ncPUU=zeta22;
nhPUU=zeta21*ncPUU;
ncIUU=zeta24;
nhIUU=zeta23*ncIUU;
lhUU=zeta20*zeta25/(zeta2+zeta20*zeta25);
lcUU=1-lhUU;
zeta26=(1-nuUU)*(1-betaGUU*(1-deltaUU))/nuUU;
zeta27=(tauhUU*zeta20+tauUU+tauUU*deltahUU*zeta20)/(zeta26+deltaUU);
zeta28=zeta1^muckUU*lcUU^muclUU*zeta27^mugUU*(ncPUU^alphaUU*ncIUU^(1-alphaUU))^(1-muckUU-muclUU-mugUU);
ycUU=zeta28^(1/(1-muckUU-mugUU));
ksUU=1/deltasUU*(zeta3*zeta20*ycUU)^muhkUU*(nhPUU^alphaUU*nhIUU^(1-alphaUU))^(1-muhkUU);
hUU=ksUU^gammaUU*lhUU^(1-gammaUU);
qlUU=zeta2*ycUU/lcUU;
qhUU=zeta20*ycUU/hUU;
qrUU=zeta12*qhUU;
cPUU=zeta18*ycUU;
lamPUU=zetaP/cPUU;
cIUU=zeta19*ycUU;
lamIUU=zetaI/cIUU;
thIUU=(1-betaIUU*RUU)*lamIUU/RUU;
cEUU=(1/xUU-zeta8-zeta9-(RUU-1)*zeta6-deltaUU*zeta1)*ycUU;
lamEUU=zetaE/cEUU;
thEUU=(1-betaEUU*RUU)*lamEUU/RUU;
wcPUU=zeta8*ycUU/ncPUU;
wcIUU=zeta9*ycUU/ncIUU;
kcUU=zeta1*ycUU;
icUU=deltaUU*kcUU;
cHUU=(deltahUU-zeta10-zeta11-(RUU-1)*zeta7-deltaUU*zeta3)*qhUU*hUU;
lamHUU=zetaH/cHUU;
thHUU=(1-betaHUU*RUU)*lamHUU/RUU;
whPUU=zeta10*qhUU*hUU/nhPUU;
whIUU=zeta11*qhUU*hUU/nhIUU;
khUU=zeta3*qhUU*hUU;
ihUU=deltaUU*khUU;
FUU=(1-1/xUU)*ycUU;
yhUU=deltahUU*hUU;
GDPUU=ycUU+qhUU*yhUU;
taxUU=tauUU*GDPUU;
kgUU=zeta27*ycUU;
gcUU=zeta26*kgUU;
gfUU=deltaUU*kgUU;
lamGUU=(1-nuUU)/gcUU;
PHIUU=lamPUU*ycUU/(1-thetaUU*betaPUU);
OMEGAUU=PHIUU;
chiHUU=gammaUU*(1-betaHUU*(1-deltahUU))/(1-betaHUU*(1-deltasUU))*lamHUU*qhUU*hUU/ksUU;
hPUU=zeta4*cPUU/qhUU;
hIUU=zeta5*cIUU/qhUU;
hRUU=zetaR*hIUU;
bIUU=mIUU/RUU*qhUU*hIUU;
bEUU=zeta6*ycUU;
bHUU=zeta7*qhUU*hUU;
dUU=bIUU+bEUU+bHUU;
bUU=dUU;
cUU=cPUU+cIUU+cEUU+cHUU;
iUU=icUU+ihUU;

lamPU =log(lamPUU);
cPU =log(cPUU);
lamIU =log(lamIUU);
cIU =log(cIUU);
thIU=log(thIUU);
bIU=log(bIUU);
hPU =log(hPUU);
hIU =log(hIUU);
hRU =log(hRUU);
qrU =log(qrUU);
qhU =log(qhUU);
wcPU =log(wcPUU);
ncPU =log(ncPUU);
wcIU =log(wcIUU);
ncIU =log(ncIUU);
whPU =log(whPUU);
nhPU =log(nhPUU);
whIU =log(whIUU);
nhIU =log(nhIUU);
RU =log(RUU);
dU =log(dUU);
FU =log(FUU);                                 
lamEU =log(lamEUU);
lamHU =log(lamHUU);
lamGU =log(lamGUU);
chiHU=log(chiHUU);
cEU =log(cEUU);
cHU=log(cHUU);
thEU =log(thEUU);
ycU =log(ycUU);
yhU=log(yhUU);
xU =log(xUU);
kcU =log(kcUU); 
thHU =log(thHUU);
khU =log(khUU);
ksU=log(ksUU);
bEU =log(bEUU);
bHU =log(bHUU);            
icU =log(icUU);
ihU =log(ihUU);                                                                                                         
gcU =log(gcUU);  
gfU =log(gfUU);  
kgU=log(kgUU);
cU =log(cUU);
iU =log(iUU);
GDPU =log(GDPUU);
bU =log(bUU);  
piU =log(piUU);
lcU =log(lcUU);
lhU =log(lhUU);
AcU =log(AcUU);
AhU =log(AhUU);
qckU =log(qckUU);
qhkU =log(qhkUU);
qlU=log(qlUU);
epscU=log(epscUU);
epshU=log(epshUU);
epslU=log(epslUU);
PHIU=log(PHIUU);
OMEGAU=log(OMEGAUU);
taxU=log(taxUU);
hU=log(hUU);
piU_ob=log(piUU);
cU_ob=0;
iU_ob=0;
GDPU_ob=0;
qhU_ob=0;
qlU_ob=0;
RU_ob=log(RUU);
bU_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
