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


%% Enter model equations here

options=optimset(); % set options for numerical solver

% the steady state computation follows FV (2006), section 4.1
PI=PIbar;
u=1;
q=1;
d=1;
df=1;
sh=1;
si=1;
phi=1;
phif=1;
m=0;
zeta=1;
mu_z=exp(LambdaYd);
mu_I=exp(Lambdamu);
mu_A=exp(LambdaA);

%set the parameter Lambdax
Lambdax=mu_z;

%set the parameter gammma1
gammma1=mu_z*mu_I/betta-(1-delta);
if gammma1<0 % parameter violates restriction; Preventing this cannot be implemented via prior restriction as it is a composite of different parameters and the valid prior region has unknown form
    check=1; %set failure indicator
    return; %return without updating steady states
end


r=1*gammma1;

R=1+(PI*mu_z/betta-1);

%set Rbar
Rbar=R;

PIstar=((1-thetap*PI^(-(1-epsilon)*(1-chi)))/(1-thetap))^(1/(1-epsilon));

mc=(epsilon-1)/epsilon*(1-betta*thetap*PI^((1-chi)*epsilon))/(1-betta*thetap*PI^(-(1-epsilon)*(1-chi)))*PIstar;

w=(1-alppha)*(mc*(alppha/r)^alppha)^(1/(1-alppha));

vp=(1-thetap)/(1-thetap*PI^((1-chi)*epsilon))*PIstar^(-epsilon);
tempvaromega=alppha/(1-alppha)*w/r*mu_z*mu_I;

% 3 Gleichungen aufstellen und mit fsolve lösen

z0=[1,1,1,1,1]
[z,fval,exitflag]=fsolve(@lddet,z0);

if exitflag <1
    %indicate the SS computation was not sucessful; this would also be detected by Dynare
    %setting the indicator here shows how to use this functionality to
    %filter out parameter draws
    check=1; %set failure indicator
    return; %return without updating steady states
end

ld_m= z(1);
ld_f = z(2);
ld_h = z(3);
wm   = z(4);
wf   = z(5);

ld = ld_h*(sh * si);

k=tempvaromega*ld;
x=(1-(1-delta)*(mu_z*mu_I)^(-1))*k;
x_h = sh^(-1) * x;
yd=(mu_A/mu_z*k^alppha*ld_h^(1-alppha))/vp;
yd_h = sh^(-1) * yd
c_h=(( mu_A/mu_z * tempvaromega^alppha)/vp - tempvaromega*(1-(1-delta)*(mu_z*mu_I)^(-1)))*ld_h*(sh * si);
lambda=1/c_h;
F=yd-1/(1-alppha)*w*ld;


g1=lambda*mc*yd/(1-betta*thetap*PI^((1-chi)*epsilon));
g2=epsilon/(epsilon-1)*g1;

c = sh * c_h;

%% end own model equations

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
