function [ys,params,check] = DJ_steadystate(ys,exo_ss,M_,options_);
% function [ys,params,check] = DJ_steadystate(ys,exo,M_,options_)
% computes the steady state for the DJ_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
%   - 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) 2013-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
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
    paramname = M_.param_names{ii};
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;
%% 
beta_t = 0.99;
alpha =1/3;
beta=0.45;
theta = 0.11;

%% 
x0 =[8;0.5];
[xx, fval, exitflag] = fsolve(@DJ_ss, x0, optimoptions('fsolve','TolFun',1e-10, 'MaxIter', 1e10),sigma_c,sigma_l,beta_t,beta,delta_k,delta_d,theta_p,alpha,theta,z,eta,epsilon_c, rhod,  rhog,rhoh,rhos);
k= xx(1);
N= xx(2);
disp('Current guess values for steady state:');
disp(xx); 
disp('Residuals of the steady state equations:');
disp(fval); 

%% Enter model equations here
d=1;
h = 1; 
g = 1; 
s =1; 
pi=1;
i=1/beta_t-1;
r=i+1;
phi_c=epsilon_c;
%% 
I= delta_k*k;
pistar=((pi^(1-eta)-theta_p)/(1-theta_p))^1/(1-eta);
mc =(1-theta_p*beta_t*pi^eta)/(1-theta_p*beta_t*pi^(eta-1))*pi/pistar*(eta-1)/eta;
R =1/beta_t-(1- delta_k);
dp= (1-theta_p)*(pi/pistar)^eta/(1- pi^eta*theta_p);
p_d=(h+g+(theta-alpha-beta)*s)/(1-alpha-beta+theta);
D=((p_d-s)/((1-alpha-beta)*d*mc*k^alpha*N^beta))^(1/(theta-alpha-beta));
A=d*D^theta;
Y= A*k^alpha * N^beta * D^(1-alpha-beta)/dp;
c=Y-I;
w=N/(phi_c*c^(-sigma_c));
x1 = c^(-sigma_c)*phi_c*mc*Y/(1-beta_t*theta_p*pi^eta);
x2=c^(-sigma_c)*phi_c*Y/(1-beta_t*theta*pi^(eta-1));


%% 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 ii = 1:NumberOfEndogenousVariables
    varname = M_.endo_names{ii};
    eval(['ys(' int2str(ii) ') = ' varname ';']);
end
