function [ys,params,check] = ghh_rbc_steadystate(ys,exo,M_,options_)
% Inputs: 
%   - ys        [vector] vector of initial values for thesteady 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] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright (C) 2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

% read out parameters to access them with their name
%BETA=NaN; %make parameter known to Matlab function, prevents crashes due to Matlab function with same name;
          %will be overwritten next
          
NumberOfParameters = M_.param_nbr;
for jj = 1:NumberOfParameters
  paramname = M_.param_names{jj};
  eval([ paramname ' = M_.params(' int2str(jj) ');']);
end

% read in instrument values
for jj = 1:size(options_.instruments,1)
  eval([options_.instruments{jj} ' = ys(strmatch(options_.instruments{jj},M_.endo_names,''exact'')) ;']);
end
% initialize indicator
check = 0;

%Fill in here ****

yoverk = (1/BETA - (1-DELTA))/ALPHA;
xoverk = DELTA;
covery = 1 - DELTA*1/yoverk-GY;
xovery = xoverk/yoverk;
kovern = (yoverk)^(1/(ALPHA-1));
phi0 = (1-ALPHA)/PHI*kovern^ALPHA*(1/3)^(1-PHI);


n = 1/3;
L = 1-n;
a = 1;
gamma = 1;
k = kovern*n;
y = a*k^ALPHA*n^(1-ALPHA);
c = covery*y;
%mu = 1/c;
mu = 1/(c - gamma*phi0*n^PHI);

w = (1-ALPHA)*y/n;
gs = 1;
yss = y;
G = gs*GY*yss;




%% end own model equations

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

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

