function [ys,params,check] = SuddenStopSGU_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

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


%% Enter model equations here

% P_N normalized to 1 
P_N = 1; 
% Eq10
P=(omegaa^(1/(1+etaa))+(1-omegaa)^(1/(1+etaa))*P_N^(etaa/(1+etaa)))^((1+etaa)/etaa);

% Solving for C_N: 

% eq7 = eq9:
% eq4 in eq11: ppi = Y_T + ((1-aplhaa)/aplhaa)*N*W and simplify. Therefore:
% C_N = (1/aplhaa)*W*N-T_N 
% -------------------------------------------------
% From eq1: N = (deltaa*(-(1/R)^(-1/B)+1+C))^(1/deltaa)
% Subs C (eq14) then C_T (eq2) in eq1: 
% eq2  C_T = ((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N
% Eq1:  N = (deltaa*(-(1/R)^(-1/B)+1+((omegaa*(((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa))))^(1/deltaa)
% -------------------------------------------------
% Subs for N in eq3:
% W= (N^(deltaa-1))*P*(1+tauu)
% W= (((deltaa*(-(1/R)^(-1/B)+1+((omegaa*(((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa))))^(1/deltaa))^(deltaa-1))*P*(1+tauu)
% -------------------------------------------------
%  -C_N+(1/aplhaa)*((((deltaa*(-(1/R)^(-1/B)+1+((omegaa*(((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa))))^(1/deltaa))^(deltaa-1))*P*(1+tauu))*((deltaa*(-(1/R)^(-1/B)+1+((omegaa*(((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa))))^(1/deltaa))-T_N 

options = optimoptions('fsolve','Display','iter','MaxIterations',400,'TolX',1e-10,'TolFun',1e-10);
[C_N,fval,exitflag]=fsolve(@(C_N)-C_N+(1/aplhaa)*((((deltaa*(-(1/R)^(-1/B)+1+((omegaa*(((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa))))^(1/deltaa))^(deltaa-1))*P*(1+tauu))*((deltaa*(-(1/R)^(-1/B)+1+((omegaa*(((1-omegaa)/omegaa)^(-1/(1+etaa))*C_N)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa))))^(1/deltaa))-T_N ...
    ,0,options);
 
if exitflag <1
    check=1;    %Set failure indicator
    disp('fsolve Failed: SS computation was not Sucessful')
    return;     %return without updating steady states
end

% if isreal(C_N)
%     disp('C_N is a real number or array.')
% else
%     disp('C_N contains complex numbers.');
% end

%Similarly: [N,fval,exitflag]=fzero(@(N) - N^(deltaa-aplhaa)+ A*aplhaa*K^(1-aplhaa)/(1+tauu),1);

% From eq2
C_T=((omegaa/(1-omegaa))^(1/(etaa+1))*C_N);
% From eq14
C = ((omegaa*(C_T)^(-etaa)+(1-omegaa)*(C_N)^(-etaa))^(-1/etaa));
% From eq1
N = (deltaa*(-(1/R)^(-1/B)+1+C))^(1/deltaa);
% From eq3:
W= (N^(deltaa-1))*P*(1+tauu); 
%From eq4
A=(P_N*aplhaa*K^(1-aplhaa)*N^(aplhaa-1)/W)^(-1);
% From eq11
ppi = Y_T+P_N*(1-aplhaa)*A*K^(1-aplhaa)*(N^aplhaa);
% From eq12
GDP = W*N + ppi;
% Use Budget Constraint eq7 for D: 
D=(W*N+ppi-(1+tauu)*(C_T + P_N*C_N)-T-P_N*T_N)/(1-R);
% Use eq13 for GT: 
GT = tauu*(C_T + P_N*C_N) + T;
% Use eq16 for u: 
u = ( (C-(N^deltaa)/deltaa)^(1-Sigmaa) -1 )/(1-Sigmaa);
% From Eq6: 
eta_m = -u/(1-(1+C-(N^deltaa)/deltaa)^(-B));
% Usin eq 15:
Uc = (C-(N^deltaa)/deltaa)^(-Sigmaa)-eta_m*(-B*(1+C-(N^deltaa)/deltaa)^(-B-1));
% Use eq5 for lamda:
lambda = Uc/(P*(1+tauu));
% Use Eq8: 
mu = 0;
% Use eq1-19:
Z=1;Z2=1;Z3=1;

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


