function [ys,params,check] = Yashiv_2022_noexp_test_steadystate(ys,exo,M_,options_)
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] 0 if steady state computation worked and to
%                        1 of not (allows to impose restrictions on parameters)

%% Step 0: initialize indicator and set options for numerical solver
check = 0;
options = optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

%% Step 1: read out parameters to access them with their name
for ii = 1:M_.param_nbr
  eval([ M_.param_names{ii} ' = M_.params(' int2str(ii) ');']);
end

%% Step 2: Check parameter restrictions
if betta<1  % parameter violates restriction (here it is artifical)
    check=1; %set failure indicator
    params=M_.params;
    return;  %return without updating steady states
end
%% Step 3: Enter model equations here
  nt = 0.889065363508680;
  xt = 0.502439035179926;
  kt = 24.613012564644315;
  ct = 2.124762549362781;
  pit = 1.000000000000000;
  at = 1;
  ut = 1-nt;
  it = 1/betta;
  %xt = (1-xt)*ut*deltan*nt;
  ht =  deltan*nt; 
  ft = at*(nt^alppha)*kt^(1-alppha); 
  gt = (e/2*(deltan)^2)*ft;
  netyt = ft+gt;
  dght =  ft/nt*e*(ht/nt);
  dgnt = -e*ft/nt*(ht/nt)^2 + alppha*gt/nt;
  dgkt = (1-alppha)/kt*gt;
  fkt = (1-alppha)*ft/kt;
  fnt = alppha*ft/nt;
  invt = deltak*kt;
  %kt = (1-alppha)*ft/((it-1+deltak)/mct+dgkt);
  %nt = alppha*ft/(dght+dgnt+wrt/mct-(1-deltan)/it*dght);
  %netyt = ct+invt;
  mct = (epsilon-1)/epsilon;
  wrt = gammma*mct*(fnt-dgnt)+(1-gammma)*(chi*ct*nt^phi+xt/(1-xt)*gammma/(1-gammma)*(mct*dght));
  mpt = 1;
  bkt = (fkt-dgkt)*mct;

%% Step 4: Update parameters and variables
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

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