function [ys,check] = Version4_28Feb_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
a=0;
sbyn=(log(phi/chi))/neta;
dn=log((chi*exp(neta*sbyn))+(1-phi));
m=log(disc)+(dn*((1-gammat)/theta));
newk = log(exp(dn)+delta-1);
newkprime=log(1);
k_h1t=log(0.6);
k_h2t=log(0.4);
bt=log(0.95);
[theta_t,fval,exitflag]=fsolve(@(theta_t)(exp(m+dn)*((1-s)+(((exp(k_h1t)*(1+exp(iota*theta_t))^(1/iota))+exp(k_h2t))^(-1))*(((1-eta)*(1-epsitlon)*(1-alphat)*exp(alphat*(log((delta-1+exp(-m))/(((epsitlon*nu)^(epsitlon/(1-epsitlon)))*(1-epsitlon)*alphat*exp((a+log((1+s*(1+exp(-iota*theta_t))^(-1/iota))))*(1-alphat))))/(alphat-1)-log((1+s*(1+exp(-iota*theta_t))^(-1/iota)))))*exp((1-alphat)*a)*((epsitlon*nu)^(epsitlon/(1-epsitlon))))-eta*(exp(k_h1t+theta_t)+ (exp(k_h2t)/((1+exp(-iota*theta_t))^(1/iota))))-(1-eta)*exp(bt))))-1,log(2),options);
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

lab=log((1+s*(1+exp(-iota*theta_t))^(-1/iota)));
ka=log((delta-1+exp(-m))/(((epsitlon*nu)^(epsitlon/(1-epsitlon)))*(1-epsitlon)*alphat*exp((a+log((1+s*(1+exp(-iota*theta_t))^(-1/iota))))*(1-alphat))))/(alphat-1);
que=log((1 + exp(theta_t*iota))^(-1/iota));
V=theta_t+log(1-exp(lab));
k_ht=log(exp(k_h1t)+exp(k_h2t+que));
ya=((alphat*ka)+((a+lab)*(1-alphat)))+(log(epsitlon*nu)*(epsitlon/(1-epsitlon)));
x=(alphat*ka)+((a+lab)*(1-alphat))+(log(epsitlon*nu)*(1/(1-epsitlon)));
pi=x+log((1/nu)-1);
firm_val=log(exp(pi)+(((1-phi)/chi)*exp((1-neta)*sbyn)));
ia   = log( exp(ka)*(exp(dn)+delta-1) );
wa=log(eta*((1-epsitlon)*(1-alphat)*exp(ya-lab) + exp(k_ht+theta))+(1-eta)*exp(bt));
ca=log(exp(ya)-exp(ia)-exp(x)-exp(sbyn)-exp(k_ht+V));
Daa=log(exp(ca)-exp(wa+lab));
Dfirm=log(exp(Daa)-exp(pi)+exp(sbyn));
rf=-m;
uc   = log(((1-disc)/(1-disc*exp((1-1/psit)*dn)))^(psit/(psit-1)));
Q    = exp((uc+dn)*(1-gammat));
logQ = log(Q);
dc=dn;

%% 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