function [y_ss Out_] = steadystate_rbc_gss_tr(S_)


%% Define parameters

alpha     = S_.param.alpha;
z_ss      = S_.param.z_ss ;
A_ss      = S_.param.A_ss ;
beta      = S_.param.beta;
delta     = S_.param.delta;
gamma     = S_.param.gamma;
nu        = S_.param.nu;
h         = S_.param.h;
by_F_ss   = S_.param.by_F_ss;
zeta_b_ss = S_.param.zeta_b_ss;
by_ss     = S_.param.by_ss;
g_ss      = S_.param.g_ss;
tau_c_ss  = S_.param.tau_c_ss;
tau_w_ss  = S_.param.tau_w_ss;
tau_K_ss  = S_.param.tau_K_ss;
sigma     = S_.param.sigma;
kappa_b   = S_.param.kappa_b;

%% Initial values

x1 = S_.initvals ;

%% Call function
options = optimset;
options = optimset(options,'MaxFunEvals',1000,'MaxIter',200,'TolFun',1e-6,'Display','Off');

F = getvals(x1);

if sum(isnan(F))>0 || sum(isinf(F))>0
    disp 'Steady state does not exist to specified iterations.'
    y_ss = [];
    Out_.iterations = [];
    return
end
  
[y_ss,fval,exitflag,output] = fsolve(@getvals,x1,options);

F = getvals(y_ss);

Out_.iterations = output.iterations;
Out_.fval       = fval; 
Out_.exitflag   = exitflag; 
Out_.F          = F ; 

%

function F = getvals(x1); 

% Variables
x = exp(x1); 

j_ = 1; 
c      = x(j_);  j_ = j_ + 1;
invest = x(j_);  j_ = j_ + 1;
lambda = x(j_);  j_ = j_ + 1;
r_F    = x(j_);  j_ = j_ + 1;
kk     = x(j_);  j_ = j_ + 1;
l      = x(j_);  j_ = j_ + 1;
w      = x(j_);  j_ = j_ + 1;
r_K    = x(j_);  j_ = j_ + 1;
gdp    = x(j_);  j_ = j_ + 1;
nx     = x1(j_); j_ = j_ + 1;  % nx can be negative
ca     = x1(j_); j_ = j_ + 1;  % ca can be negative
r_star = x(j_);  j_ = j_ + 1;
b_F    = x1(j_); j_ = j_ + 1;  % b_star can be negative
r      = x(j_);  j_ = j_ + 1;
b      = x1(j_); j_ = j_ + 1;  % b can be negative 
gy     = x(j_);  j_ = j_ + 1;   
tr     = x1(j_); j_ = j_ + 1;  % tr can be negative
p_b    = x(j_);

%% System of steady-state equations
% NB: All equations equal zero

F = [(z_ss^sigma -h*beta)/(c^sigma *(z_ss-h)^sigma) - lambda*(1+tau_c_ss); ... 
    
     r_F -  z_ss^sigma / beta;... 
     
     r -  z_ss^sigma / beta;... 
     
     p_b - (1+kappa_b*p_b)/r;...

     beta/(z_ss^sigma)*(1-delta) - 1 + beta/(z_ss^sigma)*r_K*(1-tau_K_ss);...
     
     gamma*l^nu - lambda*w*(1-tau_w_ss);...
     
     kk - (1-delta)*kk/z_ss - invest ;...
     
     A_ss*kk^(alpha)*l^(1-alpha)*z_ss^(-alpha)  - gdp ; ...
     
     r_K - alpha*gdp*z_ss/kk;...
     
     w - (1-alpha)*gdp/l;...
     
     nx - gdp + c + invest + g_ss;...
     
     ca - nx - (r_F-1)*b_F/z_ss;...
     
     r_F - r_star*exp(zeta_b_ss);...
     
     b_F - r_F*b_F/z_ss - nx;...
     
     g_ss + ((1+kappa_b*p_b)/z_ss -p_b)*b - tau_c_ss*c - tau_w_ss*w*l - tau_K_ss*r_K*kk/z_ss- tr;...
     
     g_ss/gdp - gy;...
     
     b*p_b/gdp - by_ss;...
     
     b_F/gdp - by_F_ss;...
     
               
    ];


end



end