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

% Step 0: read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;                        % Number of deep parameters.
for ii = 1:NumberOfParameters                             % Loop...
  paramname = 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 th

check = 0;

q_star = 1;
Y_star = 1;
% eta_star = 1;
% M_star = 1;
Pi_star = Pi_star_bar;
PiF_star = Pi_star_bar;
LAM_star = bbeta_star;
SX_star = 0;
X_star = 0;
R_star = Pi_star/bbeta_star;
r_star = R_star/Pi_star;
rk_star = (r_star + ddelta_star -1)/(1-tau_k_star);
mc_star = (varepsilon_star -1)/varepsilon_star + (1-bbeta_star)/varepsilon_star*gamma_p_star*(PiF_star -1)*PiF_star;
K_star = (aalpha_star*mc_star*Y_star)/rk_star;
N_star = 1/3;
Nr_star = (1/(zzeta_star*(1-chi_star)))^(1/(1+varphi_s_star));
w_star = ((1-aalpha_star)*Y_star*mc_star)/N_star;
Cr_star = w_star*(1-tau_w_star)*Nr_star;
Ns_star = N_star/(1-lambda_star) - (lambda_star)/(1-lambda_star)*Nr_star;
I_star = ddelta_star*K_star;
Gamma_star = (gamma_p_star/2)*(PiF_star -1)^2*Y_star;
d_star = (1-mc_star)*Y_star - Gamma_star;
d_s_star = lambda_star/(1-lambda_star)*d_star;
BGF_star = bg_star*Y_star;
BGF_s_star = lambda_star/(1-lambda_star)*BGF_star;
Cs_star = w_star*(1-tau_w_star)*Ns_star + (R_star/Pi_star -1)*BGF_s_star + (1-tau_d_star)*d_s_star;
C_star = lambda_star*Cr_star + (1-lambda_star)*Cs_star;
G_star = Y_star - C_star - I_star - Gamma_star;
J_star = G_star - w_star*N_star*tau_w_star - rk_star*K_star*tau_k_star -tau_d_star*d_star;
Z_star = (Y_star/(K_star^aalpha_star*N_star^(1-aalpha_star)))^(1/(1-aalpha_star));
% Us_star = (((1-chi_star)*Cs_star)^(1-sigma_s_star))/(1-sigma_s_star)*exp(zzeta_star*((sigma_s_star -1)/(1+varphi_s_star))*Ns_star^(1+varphi_s_star));
% Ur_star = (((1-chi_star)*Cr_star)^(1-sigma_r_star))/(1-sigma_r_star)*exp(zzeta_star*((sigma_r_star -1)/(1+varphi_s_star))*Nr_star^(1+varphi_r_star));
% U_star = lambda_star*Ur_star + (1-lambda_star)*Us_star;
UCs_star = ((1-chi_star)*Cs_star)^(-sigma_s_star)*exp(zzeta_star*((sigma_s_star -1)/(1+varphi_s_star))*Ns_star^(1+varphi_s_star));
% Vs_star = Us_star + bbeta_star*Us_star;
% Vr_star = Ur_star + bbeta_star*Ur_star;
% V_star = lambda_star*Vr_star + (1-lambda_star)*Vs_star;

% Rescaled variables
% Ylog_star = log(Y_star);          
% Clog_star = log(C_star);
% Pilog_star = log(Pi_star);
% Ilog_star = log(I_star);
% wlog_star = log(w_star);
% Nlog_star = log(N_star);
% Nrlog_star = log(Nr_star);
% Nslog_star = log(Ns_star);

%% Steady-state values (for policy part)
R_star_ss=R_star;
Pi_star_ss=Pi_star;
Y_star_ss=Y_star;
% C_star_ss=C_star;
% I_star_ss=I_star;
G_star_ss=G_star;
Z_star_ss=Z_star;  
BGF_star_ss=BGF_star;   
% w_star_ss =w_star;
% N_star_ss=N_star;
% Nr_star_ss=Nr_star;
% Ns_star_ss=Ns_star;


%% DO NOT CHANGE THIS PART.
%% Step 4: Update parameters and variables
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