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

% 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

options=optimset(); % set options for numerical solver

%% Enter model equations here
% Enter steady state sequentially
n_t = 1;

S_TC_ss = 1;
S_TCc_ss = 1;
S_TCI_ss = 1;

S_TC_t = S_TC_ss;
S_TCc_t = S_TCc_ss;
S_TCI_t = S_TCI_ss;

i_t = i_ss;
a_ee = 1; 
a_t = a_ee;
pi_ee = betta*(1+i_t)-1;
pi_t = pi_ee;
pi_ast_ss = (((1+pi_t)^(1-eps_p)-phi_p)/(1-phi_p))^(1/(1-eps_p)) - 1;
pi_ast_t = pi_ast_ss;

% Invesment to GDP ratio in steady state
I_Y = 0.201766;
% Liquidity demand coefficient in steady (auxiliar)
A_bar = S_TC_ss*S_TCc_ss*Ac;
% "Aggregate" money demand (auxiliar)
v = sqrt((2*Bc+i_t/(1+i_t))/(A_bar));
vI_0 = 0.01; % vI_t initval
vI_t = vI_steady(vI_0, I_Y, Ac, Bc, v, A_bar); % Solve numerically for vI_t
vc_t = sqrt(v^2-vI_t^2);
sc_t = S_TC_t*S_TCc_t*Ac*vc_t + Bc/vc_t - 2*sqrt(Ac*Bc);
sI_t = S_TC_t*S_TCI_t*AI*vI_t + BI/vI_t - 2*sqrt(AI*BI);
mc_t = (1+pi_ast_t)/(1+pi_t)*(eps_p-1)/(eps_p)*(1-betta*phi_p*(1+pi_t)^(eps_p))/(1-betta*phi_p*(1+pi_t)^(eps_p-1));
r_t = ((1+2*(S_TC_t*S_TCI_t*AI*vI_t-sqrt(Ac*Bc)))*(1-betta*(1-delta)))/betta;
ups_t = (1-phi_p)/(1-phi_p*(1+pi_t)^eps_p)*((1+pi_t)/(1+pi_ast_t))^eps_p;
k_t = ((1-alpa)*mc_t*a_t*(n_t)^(alpa)/(ups_t*r_t))^(1/alpa);

I_t = delta*k_t;
y_ss = (I_Y/I_t)^(-1);
y_t = y_ss;
alpa = 1-k_t*r_t/(mc_t*y_t);
c_t = (y_t-I_t*(1+sI_t))/(1+sc_t);
lambda_t = c_t^(-sig)/(1+2*(A_bar*vc_t-sqrt(Ac*Bc)));
m_t = c_t/vc_t;
w_t = (alpa)*mc_t*y_t/n_t;

mu_t = lambda_t*(1+2*(A_bar*vI_t-sqrt(AI*BI)));
theta = w_t*lambda_t/n_t^(xi);

xn_t = c_t^(-sig)*y_t*mc_t/(1-betta*phi_p*(1+pi_t)^eps_p);
xd_t = c_t^(-sig)*y_t/(1-betta*phi_p*(1+pi_t)^(eps_p-1));


% check actual values
% ups_t
% y_t
% alpa
% a_t*(k_t^(1-alpa))*(n_t^alpa)/ups_t
% phi_p
% pi_t
% phi_p*(1+pi_t)^eps_p
% mc_t

% residual for a marginal cost equation would be:
w_t*n_t + r_t*k_t - mc_t*y_t

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