%EXOG       plot results from exog.f.
%

%           Ellen McGrattan, 4-22-98
%           Revised, 8-17-98, ERM


%
% Load data.
%

load exog.dat
X     = exog;
clear exog

%
% Assign parameters and the hours function. 
%

a     = X(1);
alpha = X(2);
beta  = X(3);
delta = X(4);
gamma = X(5);
glab  = X(6);
gtec  = X(7);
kbar  = X(8);
psi   = X(9);
sigma = X(10);
ptr   = 11;
nx    = X(ptr);
ptr   = ptr+1;
xa    = X(ptr:ptr+nx-1);
ptr   = ptr+nx;
ns    = X(ptr);
ptr   = ptr+1;
taux  = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
nnzp  = X(ptr);
ptr   = ptr+1;
ia    = X(ptr:ptr+nnzp-1);
ptr   = ptr+nnzp;
ja    = X(ptr:ptr+nnzp-1);
ptr   = ptr+nnzp;
pi    = sparse(ia,ja,X(ptr:ptr+nnzp-1));
ptr   = ptr+nnzp;
n     = nx*ns;
c     = zeros(nx,ns);
c(:)  = X(ptr:ptr+n-1);
c     = c*kbar;
gama  = gamma^(alpha-1);
gam1  = (1-gamma^alpha)/(1-gamma);
alph1 = 1-alpha;
delt1 = 1-delta;
sigm1 = 1-sigma;
betg  = beta*(1+gtec)^(-sigma);
gtec1 = gtec+1;
glab1 = glab+1;
x     = xa*ones(1,ns);


%
% Construct other decision functions.
%


plotdf = input('Plot decision functions? (1=yes, 0=no) ');

if plotdf;
  ka    = kbar*xa;
  taux1 = ones(nx,1)*taux';
  k     = ka*ones(1,ns);
  h     = zeros(nx,ns);

  for j=1:ns;
    for i=1:nx;
      c0   = c(i,j);
      k0   = ka(i);
      tem1 = (1-alpha)*a*k0^alpha;
      tem2 = tem1/psi/c0;
      h0   = tem2/(tem2+gama);
      if (h0>gamma);
        h0 = (tem2-(gama-gam1)*h0)/(tem2+gam1);
      end;
      lei  = 1-h0;
      for l=1:20
        tem   = psi*c0 - tem1*h0^(-alpha)*(1-h0);
        dtem  = tem1*h0^(-alpha)*(alpha*(1-h0)/h0+1);
        h0    = h0 - tem/dtem;
        lei   = 1-h0;
      end;
      h(i,j) = h0;
    end;
  end;
  y      = a*k.^alpha.*h.^(1-alpha);
  kp     = (y-c+delt1*k)/glab1/gtec1;
  xk     = glab1*gtec1*kp-delt1*k;

  plot(xa,c),grid, title('consumption'),pause
  plot(xa,h),grid, title('hours of work'),pause
  plot(xa,xk),grid, title('investment in physical capital'),pause
  plot(xa,kp/kbar,xa,xa),grid, title('capital next period'),pause
end;

%
% Compute steady states.
%

compss = input('Compute steady state? (1=yes, 0=no) ');
if compss;
  istate  = input('Which state for policies? ');
  tau     = taux(istate);
  hss     = .1;
  for i=1:20;
    kss   = ( (1+tau)*(1-betg*(1-delta)) / (betg*alpha*a) )^(1/(alpha-1))*hss;
    xss   = (gtec1*glab1-1+delta) * kss;
    yss   = a*kss^alpha*hss^(1-alpha);
    css   = yss - xss;
    tem   = psi*css - a*(1-alpha)*(kss/hss)^alpha*(1-hss);
    dtem  = a*(1-alpha)*(kss/hss)^alpha*(alpha*(1-hss)/hss+1);
    hss   = hss - tem/dtem;
  end 
  disp(' ')
  disp(' ')
  disp(' ')
  disp('Steady state values for Exogenous Growth Model')
  disp('----------------------------------------------')
  disp(sprintf('   h                 =   %7.4f ',hss))
  disp(sprintf('   c/y               =   %7.4f ',css/yss))
  disp(sprintf('   x/y               =   %7.4f ',xss/yss))
  disp(sprintf('   k/y               =   %7.4f ',kss/yss))
  disp('----------------------------------------------')
  disp(' ')
  disp(' ')
  disp(' ')
end;

%
% Simulate time series.
%

plotts = input('Simulate and plot time series? (1=yes, 0=no) ');

if plotts;
  t0     = input('Initial date (<= 1820)?           ');
  x0     = input('Initial k[0]/kbar (in [0,1])?     ');
  i0     = input('Initial fiscal state?             ');
  nc     = input('Number of countries?              ');
  nt     = 1990-t0+1;

  Ct     = zeros(nc,nt);
  Xkt    = zeros(nc,nt);
  Gt     = zeros(nc,nt);
  Yt     = zeros(nc,nt);
  Kt     = zeros(nc,nt);
  St     = zeros(nc,nt);
  for i=1:nc;
    disp(sprintf('Computing statistics for country %g',i))
    xt(1)  = x0;
    st     = markov(full(pi),nt,i0,i);
    taux1t = 1+taux(st);
    for t=1:nt;
      i1         = find(x(:,1)<xt(t));
      if isempty(i1);  i1=1;   end;
      i1         = i1(length(i1));
      if (i1>nx-1); i1=nx-1; end;
      x1         = x(i1);
      x2         = x(i1+1);
      dx         = x2-x1;
      xi         = 2*(xt(t)-x1)/dx-1;
      bs1        = 0.5*(1-xi);
      bs2        = 0.5*(1+xi);
      ct(t,1)    = (c(i1,st(t)) *bs1 +c(i1+1,st(t)) *bs2 )/kbar;

      tem1       = alph1*a*(xt(t)*kbar)^alpha / kbar;
      tem2       = tem1/psi/ct(t);
      h0         = tem2/(tem2+gama);
      if (h0 > gamma);
        h        = (tem2-(gama-gam1)*gamma)/(tem2+gam1);
      end
      lei        = 1-h0;
      for k=1:10;
        tem      = psi*ct(t) - tem1*h0^(-alpha)*lei;
        dtem     = tem1*h0^(-alpha)*(alpha*lei/h0+1);
        h0       = h0 - tem/dtem;
        lei      = 1-h0;
      end;
      ht(t,1)    = h0;
      yt(t,1)    = a*(xt(t)*kbar)^alpha*ht(t)^alph1/kbar;
      xkt(t,1)   = yt(t)-ct(t);
      xt(t+1,1)  = (delt1*xt(t)+xkt(t))/glab1/gtec1;
    end;
    xt       = xt(1:nt);
    Ct(i,:)  = ct' * kbar;
    Lt(i,:)  = ht';
    Xkt(i,:) = xkt' * kbar;
    Yt(i,:)  = yt' * kbar;
    Kt(i,:)  = xt' * kbar;
    St(i,:)  = st';
  end;
  At       = (ones(nc,1)*gtec1.^[1:nt]);
  Nt       = (ones(nc,1)*glab1.^[1:nt]);     % Total population
  Lt       = Lt.*Nt;                         % Total number of workers
  Yt       = Yt.*At.*Nt;
  Ct       = Ct.*At.*Nt;
  Kt       = Kt.*At.*Nt;
  Xkt      = Xkt.*At.*Nt;
  yl       = Yt./Lt;
  yn       = Yt./Nt;
  ylrel    = yl./(ones(nc,1)*exp(mean(log(yl))));
  ynrel    = yn./(ones(nc,1)*exp(mean(log(yn))));
  syn      = sort(yn);
  g1       = ((yl(:,1972-t0+1)./yl(:,1960-t0+1)).^(1/12)-1)*100;
  g2       = ((yl(:,1985-t0+1)./yl(:,1972-t0+1)).^(1/13)-1)*100;
  gy       = ((yl(:,1985-t0+1)./yl(:,1960-t0+1)).^(1/25)-1)*100;
  gy2      = 1/10*(log(yn(:,11:10:nt))-log(yn(:,1:10:nt-10)));
  yrs      = [1820,1870,1890,1913,1950,1973,1989];
  T        = t0:t0+nt-1;

  figure(1)
  semilogy(yrs,syn(250,yrs-t0+1)*1000/syn(500,1820-t0+1),'+', ...
           yrs,syn(750,yrs-t0+1)*1000/syn(500,1820-t0+1),'+', ...
           yrs,syn(100,yrs-t0+1)*1000/syn(500,1820-t0+1),'o', ...
           yrs,syn(900,yrs-t0+1)*1000/syn(500,1820-t0+1),'o')
  axis([1800,2000,2*10^2,3*10^5]);
  title('Fig. 1. GDP Per Capita, 1800-1989')
  ylabel('GDP Per Capita')

  figure(2)
  hist0(ynrel(:,1989-t0+1));
  axis([log(1/8),log(8),0,.39])
  title('1989 Relative GDP Per Capita (World Average = 1)')
  ylabel('Fraction of Countries (%)')

  figure(3)
  plot(log(ylrel(:,1960-t0+1)),gy,'o'), ...
  axis([log(1/8),log(8),-3,8])
  set(gca,'XTickLabel',['1/8';'1/4';'1/2';' 1 ';' 2 ';' 4 ';' 8 '], ...
          'XTick',log([1/8;1/4;1/2;1;2;4;8]))
  title('Growth versus Initial GDP Per Worker, 1960-1985')
  xlabel('Relative GDP Per Worker in 1960')
  ylabel('Growth Rate of GDP Per Worker, 1960-85')

  figure(4)
  plot(g1,g2,'o')
  axis([-5,10,-5,10])
  title('Persistence of Growth Rates, 1960-1985')
  xlabel('Growth rates of GDP per Worker, 1961-72')
  ylabel('Growth rates of GDP per Worker, 1973-85')
  
  figure(5)
  sgy2 = sort(gy2);
  plot(T(11:10:nt),100*mean(sgy2(976:1000,:)))
  axis([1880,2000,2,12])
  title('Maximum GDP Per Capita Growth, 1870-1990')
  ylabel(' Maximum Average Annual Growth Rate (%) ')


  %
  % Data for figures in fig1a,fig1b,fig2,fig3,fig4
  %

  fig1a = [yrs',syn(250,yrs-t0+1)'*1000/syn(500,1820-t0+1);
           yrs',syn(750,yrs-t0+1)'*1000/syn(500,1820-t0+1)];
  fig1b = [yrs',syn(100,yrs-t0+1)'*1000/syn(500,1820-t0+1);
           yrs',syn(900,yrs-t0+1)'*1000/syn(500,1820-t0+1)];
  fig2  = [log(ylrel(:,1960-t0+1)),gy];
  fig3  = [g1,g2];
  fig4  = [T(11:10:nt)',100*mean(sgy2(976:1000,:))'];

  %
  % Mobility matrix as in Chari-Kehoe-McGrattan:
  %
  mob  = zeros(6);
  bin  = [0,1/4,1/2,1,2,4,10000];
  for i=1:6;
    ii = find(ylrel(:,1960-t0+1) > bin(i) & ylrel(:,1960-t0+1)<= bin(i+1));
    for j=1:6;
      jj = (ylrel(ii,1985-t0+1) > bin(j) & ylrel(ii,1985-t0+1)<= bin(j+1));
      if ~isempty(jj);
        mob(i,j)  = sum(jj);
      end;
    end
  end;

  mob = mob./(sum(mob')'*ones(1,6));
  disp('Mobility Matrix:')
  disp(sprintf(' %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',round(mob*100)'/100))

end;

