function [ys,check] = drh1nk_steadystate(ys,exo)
% function [ys,check] = alpha_steadystate(ys,exo)
% computes the steady state for the alpha.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;

%% 
beta=0.99;
delta=0.025;
sigma=2;
eta=1.2;
h=0.6;
alphal=0.6;
alphak=1-alphak;
tau=2;
thetaw=0.75;
thetap=0.75;
phi= 1;
psi=1;
phipi = 1.5;
phiy = 0.5;
rhoak =0.9;
rhoal =0.9;
rhomw =0.9;
rhomp =0.9;
rhoi =0.9;
rhoupsilon =0.9;
rhouv =0.9;
sdg =0.01;
sdak =0.01;
sdal =0.01;
sdmw =0.01;
sdmp =0.01;
sduv =0.01; 
sdi =0.01;
sdupsilon =0.01;
chi1 = 1/beta - (1-delta);
chi2 = 0.01;
YSS     = 1.1796;
ISS     = 1/beta;
PISS    = 1;  



uv=1;
upsilon=1;
Ak=1;
Al=1;
DS  =0.1;
LS =(1-DS)*alphal;
KS =1-LS-DS;
Al=1;
u=1;
vp = 1;
pi=1;
i=pi/beta;
pistar=1;
Alc=1;
Akc=1;
mwc=1;
mpc=1;
upsilonc=1;
zak=0;
zal=0;
zmp=0;
zmw=0;
zupsilon=0;
%%
x0 = [1];
[xx,fval]=fsolve(@nkss,x0,optimoptions('fsolve','TolFun',1e-5, 'MaxIter', 1e15),eta,beta,delta,h,sigma,phi,psi,thetaw,thetap,upsilon,DS,LS,KS,alphal,alphak);
%[xx,fval]=fsolve(@nkss,x0,optimoptions('fsolve','TolFun',1e-5, 'MaxIter', 1e15));
L = xx(1);


Rk= (1/beta - (1-delta))/upsilon;
mp=alphak*KS^(-1/eta)*Rk^((1-eta)/eta);
Y=(alphal/mp/LS)^(eta/(eta-1))*L;
W=LS*Y/L;
K=KS*Y/Rk;
Inv=delta*K;
X = Inv/upsilon;
C = Y - X;
D=DS*Y;
Wstar = W;
lambda = (C-h*C)^(-sigma)-beta*h*(C-h*C)^(-sigma);
Kbar = K;
mu = lambda/upsilon;
Q  = mu/lambda;
H1 = psi*L^(1+phi)/(1-thetaw*beta);
H2 = lambda*L/(1-thetaw*beta);
mc= 1/mp;
X1 = lambda*mc*Y/(1-thetap*beta);
X2 = lambda*Y/(1-thetap*beta);
epsilonp = mp/(mp-1);
%mw = 1/W/mp;
%mw=mp;
mw=W*lambda/psi/L^phi;
epsilonw = mw/(mw-1);
Alt=Al;
Akt=Ak;
mpt=mp;
mwt=mw;
upsilont=upsilon;
%% 

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
