function [ys,check] = thesis_steadystate(ys,exe)
% 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/>.

global M_

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for ii = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(ii,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(ii) ');']);         %    Get the value of parameter i.
end                                                           % End of the loop.  
check = 0;



%% Enter model equations here


% the steady state computation 
% Sss= 1; 
% Ystarss=1; 
% Cstarss = 1; 
% NEXss = 0; 
% Pss = 1; 
% PIss = 1; 
% Pfss = 1; 



% variables that don't depend on Iss:
Phss = Pfss/Sss; 
% Pss = ((1-gammaa)*(Phss)^(1-a) + gammaa*(Pfss)^(1-a))^(1/(1-a)); % Pss = Phss = Pfss 
% PIss = ((1-gammaaI)*(Phss)^(1-aI) + gammaaI*(Pfss)^(1-aI))^(1/(1-aI)); % PIss = Phss = Pfss 
Qss = Pfss/Pss; % Qss = 1 
EXss = (Phss/Pfss)^(-aW)*Ystarss; % EXss =1
Mss = EXss - NEXss; % Mss = 1 
MCss = (tet-1)/tet; 
Rkrss = (1-(betaa*(1-deltaa)))/betaa;
Wrss = ((tet-1)/tet)^(1/(1-alphaa))*(1-alphaa)*(alphaa*betaa/(1-(betaa*(1-deltaa))))^(alphaa/(1-alphaa));


% 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 directly or indirectly: 

Nss = (((1-(gammaaI*Iss))/gammaa) + 1)/(alphaa*((tet-1)/tet)* betaa/(1-(betaa*(1-deltaa))))^(alphaa/(1-alphaa)); 
Css = (1-(gammaaI*Iss))/gammaa;
Yss = ((1-(gammaaI*Iss))/gammaa) + 1; 
Kss = (alphaa*((tet-1)/tet)* betaa/(1-(betaa*(1-deltaa))))^(1/(1-alphaa))*Nss;
Aoness = Yss*(tet-1)/tet*(1-(omegaa*betaa)); 
Atwoss= Yss/(1-(omegaa*betaa)); 
vss = (Cstarss/Css)^(-sigmaa)*1/Qss; 



%% end own model equations


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





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