% program to solve the deterministic growth model by global projection, using polynomial approximation
% NOT A NUMERICALLY EFFICIENT OR STABLE IMPLEMENTATION, just a proof of concept
% Michael Reiter, February 2020
% Input: order of the polynomial approximation
function [coef,res] = growthnonlin(n)

  delta = 0.025;
  alpha = 0.4;
  beta = 0.99;

  % solve for steady state:
  % 1 = beta*(1+alpha*k.^(alpha-1) - delta)
  rstst = 1/beta-1;
  global kstst;
  kstst = (alpha/(rstst+delta))^(1/(1-alpha));

  % grid of capital at which Euler residual will be set to zero:
  width = 0.2*kstst;  % approximate within plus/minus 20 percent of steady state
  gridK = linspace(kstst-width,kstst+width,n)';
  Invstst = delta*kstst;
  cstst = kstst^alpha - Invstst;

  % initial guess: c = cstst + 0.02*(k-kstst)
  % constant is cstst, since k is measured in deviations from steady state:
  coefStart = [zeros(n-2,1);0.02;cstst];
  res0 = eulerresid(gridK,coefStart);
  disp(res0);
  func = @(coef) eulerresid(gridK,coef);
  %   coef = fsolve(func,coefStart);
  coef = broydn(func,coefStart,[1e-8;0;1]);

  % check accuracy: compute Euler residual at many grid points
  gridKbig = linspace(kstst-width,kstst+width,1000)';
  res = eulerresid(gridKbig,coef) / util_c(cstst);
end



% coefficients of polynomial are defined for capital in deviations from steady state:
function c = consfunc(k,coef)
  global kstst;
  kdev = k-kstst;
  c = polyval(coef,kdev);
end

function muc = util_c(c)
  gam = 1;
  muc = c.^(-gam);
end

function resid = eulerresid(kbeg,coef) % k is capital at the beginning of the period
  delta = 0.025;
  alpha = 0.4;
  beta = 0.99;
  c = consfunc(kbeg,coef);
  if any(c<=0)
    resid = 1e100;  % signal to broydn that this is an inadmissible coefficient vector
    return
  end
  muc = util_c(c);
  kend = (1-delta)*kbeg + kbeg.^alpha - c;
  cnext = consfunc(kend,coef);
  if any(cnext <=0) || any(kend <=0)
    resid = 1e100;
    return
  end
  mucnext = util_c(cnext);
  rnext = alpha*kend.^(alpha-1) - delta;
  resid = beta*(1+rnext).*mucnext - muc;
end
