% computes the steady state of BNT (Blanchard-Yaari model) with growth (per capita version) by search
% Rob Luginbuhl January 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;
%  disp('gamma');
%  disp(pgamma);
%  disp('beta');
%  disp(pbeta);
%  disp('chi');
%  disp(chi);
%  disp('theta');
%  disp(theta);
%  disp('d');
%  disp(d); 
  % have you done a grid search on r yet? (no)
  rooster = 0;
  % set step length for grid search
  rstep = 0.001;
  % technology
  dA     = exp(pgamma);
  % dividend
  divai = 0.0;
  % labor
  % lab = 1.0;
  % help variable for the disturbance
  e = 1.0;
  % Non-stationary consumption must be 1 for Dynare
  C_obs = 1.0;
  % starting guess value for interest rate: (1/pbeta -1)/2 + (1/(pbeta*d)-1)/2
  %r = (1.0 + d - 2*pbeta*d)/(pbeta*d*2);
  % try this instead: starting value for r, existing simulation value:
    % get location of parameter in global
    p_name = 'r';
    i = strmatch(p_name, M_.endo_names, 'exact');
    r = oo_.steady_state(i);
    %disp(r);
    if r <= 0
        r = (1.0 + d - 2*pbeta*d)/(pbeta*d*2);	 % must not be negative
    elseif r >= 2
        r = (1.0 + d - 2*pbeta*d)/(pbeta*d*2);   % must not be greater than 1
    end
    if r <= 0
        r = rstep;	 % still need help with good starting value for r
    elseif r >= 2
        r = 0.5;   % must not be greater than 1
    end
    %disp(r);
  % need endogenous variable outside while loop to avoid matlab crash:
  % 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;
  % savings
  n = k;
  % capital rental rate
  rk = r;
  % keep old (guess) value of interest rate
  oldr = r;
  % update value of interest rate
  newr = (c - w + k*(1.0 - 1.0/dA))*dA/k;
  if (newr >= 2)
      newr = 0.1; % got to be something reasonable, also a check for Inf
  end
  if (newr <= 0)
      newr = 0.01; % don't want negative number or zero
  end
  %disp(newr);
  % start with new and old values of r each weighted by half to get new one
  wigt = 1;
  % new guess is inbetween
  r = (wigt*oldr + newr)/(wigt + 1);
  %disp(r);
  % count the number of loops
  nloop = 1;
  % loop through until converge give r = new r
  while (round(abs((newr - oldr)/r)*1.0d13) > 0.0)
      % count the loops
      nloop = nloop + 1;
      %disp(nloop);
      % keep old (guess) value of interest rate
      oldr = r;
      % update value of interest rate
      newr = newrate(r, chi, dA, d, theta, pbeta);
      %disp(newr);
      % new guess is inbetween
      r = (wigt*oldr + newr)/(wigt + 1);
      %if (abs(r) >= Inf)
      %  disp ('trouble Inf begin!');
      %  disp(nloop);
      %  disp(d);
      %  disp(pbeta);
      %  disp(pgamma);
      %  disp(theta);
      %  disp(chi);
      %  disp(rooster);
      %  disp(wigt);
        %disp(w);
        %disp(k);
        %disp(hw);
        %disp(h);
        %disp(s);
        %disp(c);
      %  disp(r);
      %end
      %disp(r);
      if (nloop > 50) | (abs(r) >= Inf)
          % look for new starting values if rooster = 0
          if (rooster < 1)
              rooster = 2; %do not do this grid search more than once
              r = 0; % initialize interest rate for grid search
              % look for starting value under 1 (really a check for Inf)
              newr = Inf;
              while (abs(newr) >= Inf) % look for first good value
                % if newr is too big (Inf) try next value in grid
                r = r + rstep;
                % update value of interest rate
                newr = newrate(r, chi, dA, d, theta, pbeta);
                % keep last value 
                oldr = r;
              end
              % how far off from convergence are we?
              diferr = oldr - newr;
              %disp (r);
              %disp(diferr);
              % assume sign of difference is positive
              difsign = 1;
              if (diferr < 0)
                  % the sign turns out to be negative
                  difsign = -1;
              end
              % do grid search until sign of error changes
              olddif = 0; %Matlab seems to want this
              while (diferr * difsign > 0) && (r < 2.0)
                  oldr = r; % keep old value of interest rate
                  olddif = diferr; % keep error associated with oldr
                  r = r + rstep; % move along grid for interest rate
                  % update value of interest rate
                  newr = newrate(r, chi, dA, d, theta, pbeta);
                  % how far off from convergence are we?
                  diferr = r - newr;
                  %disp (r);
                  %disp(diferr);
              end
              if (olddif * diferr < 0) % sign change occured, get interpolation
              % linear interpolation of new value of interest rate
                newr = oldr - olddif * (oldr - r)/(olddif - diferr);
                r = newr;
              else
                  disp('Grid search failure!');
                  disp(chi);
                  disp(pgamma);
                  disp(pbeta);
                  disp(theta);
                  disp(d);
                  disp(r);
                  disp(diferr);
                  disp(olddif);
                  r = 0.05; % random alternative try :P
                  newr = r;
                  rstep = rstep/10; % make a finer grid for search
                  rooster = 0; % allow a new grid search
              end;
              % keep old (guess) value of interest rate
              oldr = r;
              % update value of interest rate
              newr = newrate(r, chi, dA, d, theta, pbeta);
              %disp(r);
              %disp(r-newr);
              % new guess is inbetween
              r = (wigt*oldr + newr)/(wigt + 1);
          end    
      %    disp('Trouble with Steady State!');
      %    disp(nloop);
      %    disp(d);
      %    disp(pbeta);
      %    disp(theta);
      %    disp(chi);
          disp(wigt);
          nloop = 0;
          wigt = wigt*2;
      end
  end
  % get state values given steady state r
    % 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;
  % savings
  n = k;
  % capital rental rate
  rk = r;
  if (w+k+hw+h+s+c+r >= NaN)
      disp ('trouble NaN!');
      disp(nloop);
      disp(d);
      disp(pbeta);
      disp(pgamma);
      disp(theta);
      disp(chi);
      disp(rooster);
      disp(wigt);
      disp(w);
      disp(k);
      disp(hw);
      disp(h);
      disp(s);
      disp(c);
      disp(r);
  end
  if (abs(w+k+hw+h+s+c+r) >= Inf)
      disp ('trouble Inf!');
      disp(nloop);
      disp(d);
      disp(pbeta);
      disp(pgamma);
      disp(theta);
      disp(chi);
      disp(rooster);
      disp(wigt);
      disp(w);
      disp(k);
      disp(hw);
      disp(h);
      disp(s);
      disp(c);
      disp(r);
  end
  %disp(nloop);
  % convergence was with this value of r
  % r = oldr; this is really silly if there has been convergence: just keep
  % best value: average old and new
  %disp(r);
      disp(nloop);
      disp(d);
      disp(pbeta);
      disp(pgamma);
      disp(theta);
      disp(chi);
      disp(rooster);
      disp(wigt);
      disp(w);
      disp(k);
      disp(hw);
      disp(h);
      disp(s);
      disp(c);
      disp(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 newr = newrate(r, chi, dA, d, theta, pbeta)
% 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;
% update value of interest rate
newr = (c - w + k*(1.0 - 1.0/dA))*dA/k;