function [ys,params,check] = thesis_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% 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
%   - 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 <http://www.gnu.org/licenses/>.

%% 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) ');']);         %    Get the value of parameter i.
end                                                                     % End of the loop.  




%% Step 2: Enter model equations here

% values of the parameters I use here: 
betaa = 0.99; 
tet = 6; 
alphaa = 1/3; 
sigmaa = 1.47; 
etaa = 0.028; 
gammaa = 0.33; 
gammaaI = 0.18;  
deltaa = 0.05; 

% variables that don't depend on Iss:
Rkrss = (1-(betaa*(1-deltaa)))/betaa;
Mss =1; 

% find Iss : 
[Iss,fval,exitflag] =fzero(@(Iss)(((((1-(gammaaI*Iss))/gammaa) + 1)/(alphaa*(tet-1)/tet*betaa/(1-(betaa*(1-deltaa))))^(alphaa/(1-alphaa)))^etaa /...
(((tet-1)/tet)^(1/(1-alphaa))*(1-alphaa)*(alphaa*betaa/(1-(betaa*(1-deltaa))))^(alphaa/(1-alphaa)))) - (((1-(gammaaI*Iss))/gammaa)^(-sigmaa)),0.25,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


% variables that depend on Iss: 

Css = (1-(gammaaI*Iss))/gammaa;
Yss = ((1-(gammaaI*Iss))/gammaa) + 1; 



%% Step 3: 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


% ys = zeros(M_.orig_endo_nbr,1);
% for ii = 1:M_.orig_endo_nbr                                                               
%   eval(['ys(' int2str(ii) ') = ' M_.endo_names{ii} ';']);
%   ys(ii)=0;
% end 


% NumberOfEndogenousVariables = M_.endo_nbr;                    
% ys = zeros(NumberOfEndogenousVariables,1);                  
% for ii = 1:NumberOfEndogenousVariables                         
% %   varname = deblank(M_.endo_names(ii,:));                                          
% %   eval(['ys(' int2str(ii) ') = ' varname ';']);
%   ys(ii)=0;%    Get the steady state vZNue of this variable.
% end 
