
function [ys,params,check] = epusmss_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 3: Enter model equations here
r=1/betta;
w=log(w_bar);
n=log(n_bar);
cr=((w/(n^phii))^(1/sigmma))/(1-h);
a=1;
y=n/(1-covv);

std_r=(((1-rrho_stdr)*std_ssr)+std_uncr)/(1-rrho_stdr);

std_g= (((1-rrho_stdg)*std_ssg)+std_unc)/(1-rrho_stdg)

pii=((r-std_r)/(y^upssilon_y))^(1/upssilon_pii);

g=((1-rrho_g)*(rrho_y*(1-covv)*y)+std_g)/(1-rrho_g);  


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