% computes the steady state of BNT (Blanchard-Yaari model) with growth (per capita version) 
% Rob Luginbuhl September 2010
function [ys,check] = BNTnlc_nonS_PP_steadystate(ys,exe)
  global M_ oo_
  
%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end                                                           % End of the loop.  
check = 0;
%%
%% END OF THE FIRST MODEL INDEPENDENT BLOCK

%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.
%%
%  format long;
  % technology
  dA     = exp(pgamma);
  % labor
  % lab = 1.0;
  % help variable for the disturbance
  e = 1.0;
  % Non-stationary consumption must be 1 for Dynare
  C_obs = 1.0;
% get interest rate and more
[r w k hw h s c] = newrate(chi, dA, d, theta, pbeta);
  % savings
  n = k;
  % capital rental rate
  rk = r;
%%
%% END OF THE MODEL SPECIFIC BLOCK.


  %% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state values of the endogenous variables of
%% the model.
%%
NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for i = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
  eval(['ys(' int2str(i) ') = ' varname ';']);                %    Get the steady state value of this variable.
end                                                           % End of the loop.
%%
%% END OF THE SECOND MODEL INDEPENDENT BLOCK.
%
%
%subfunction to get new value of interest rate
%
function [r,w,k,hw,h,s,c] = newrate(chi, dA, d, theta, pbeta)
% interest rate function needs to be zero
f = @(r) (   (d/(1.0 + r))^(1.0 - 1.0/theta)*(d*pbeta)^(1.0/theta) - 1 +      (1 + dA*chi/(r*(1 - chi))*((1 + r)/dA - 1)) /  (1/(1 - d*dA/(1 + r)) + (1 + r)/r*chi/(1 - chi) ));
% starting guess for interest rate
rguess = 1/pbeta - 1;
% find solution for interest rate
r = fzero(f, rguess);
% wage rate
w = (r/chi)^(chi/(chi - 1.0))*(1.0 - chi);
% capital stock
k = dA*w/r/(1.0 - chi)*chi;
% permanent income
hw = w/(1.0 - d*dA/(1.0 + r));
% total wealth
h = hw + (1.0 + r)*k/dA;
% savings rate
s = (d/(1.0 + r))^(1.0 - 1.0/theta)*(d*pbeta)^(1.0/theta);
% consumption
c = (1.0 - s)*h;