%TWOSECT    plots results for TWOSECT.F which computes an equilibrium of
%           a two-sector endogenous growth model with stochastic tax rates
%           and government spending.
%           (See twosect.f for details of the model.)

%           Ellen McGrattan, 6-26-97
%           Revised, 8-17-98, ERM

% Load solution located in twosect.dat

load twosect.dat
X     = twosect;
clear twosect

a     = X(1);
alpk  = X(2);
alph  = X(3);
b     = X(4);
beta  = X(5);
chi   = X(6);
delk  = X(7);
delh  = X(8);
eta   = X(9);
lab0  = X(10);
psi   = X(11);
sigma = X(12);
thetk = X(13);
theth = X(14);
zeta  = X(15);
ptr   = 16;
nx    = X(ptr);
ptr   = ptr+1;
xa    = X(ptr:ptr+nx-1);
ptr   = ptr+nx;
ns    = X(ptr);
ptr   = ptr+1;
tauxk = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
tauxh = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
tauk1 = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
tauk2 = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
tauh1 = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
tauh2 = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
tauc  = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
z1    = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
z2    = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
z3    = 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;
c     = zeros(nx,ns);
c(:)  = X(ptr:ptr+nx*ns-1);
ptr   = ptr+nx*ns;
v     = zeros(nx,ns);
v(:)  = X(ptr:ptr+nx*ns-1);

alp   = alpk+alph;
alp1  = 1-alp;
alpi  = 1/alp;
aa    = a^alpi;
alphak= alpk*alpi;
alphah= alph*alpi;
alphai= alp1*alpi;
delk1 = 1-delk;
delh1 = 1-delh;
sigm1 = 1-sigma;
thet  = thetk+theth;
thet1 = 1-thet;
zeta1 = 1+zeta;
sigz  = sigma-zeta*sigm1;
betn  = beta*(1+eta/100)^sigz;

%
% Construct other decision functions.
%

tauk11 = ones(nx,1)*(1-tauk1)';
tauk21 = ones(nx,1)*(1-tauk2)';
tauh11 = ones(nx,1)*(1-tauh1)';
tauh21 = ones(nx,1)*(1-tauh2)';
tauc1  = ones(nx,1)*(1+tauc)';
tauxh1 = ones(nx,1)*(1+tauxh)';
tauxk1 = ones(nx,1)*(1+tauxk)';
z1     = ones(nx,1)*z1';
z2     = ones(nx,1)*z2';
z3     = ones(nx,1)*z3';
x      = xa*ones(1,ns);

for j=1:ns
  for i=1:nx
    phi   = tauk11(i,j)*tauh21(i,j)*alpk*theth/(tauk21(i,j)*tauh11(i,j)* ...
            alph*thetk);
    u(i,j)= 1/(1-phi+phi/v(i,j));
    tem1  = psi*c(i,j);
    tem2  = tauh11(i,j)*alph*aa*(x(i,j)*v(i,j)/u(i,j))^alphak*z2(i,j)^ ...
            alphai/tauc1(i,j);
    lab   = lab0;
    lei   = 1-lab;
    for k=1:10;
      tem  = tem1 - tem2*lab^(-alphak)*lei;
      dtem = tem2*lab^(-alphak)*(alphak*lei/lab+1);
      lab  = lab - tem/dtem;
      lei  = 1-lab;
    end
    l(i,j) = lab;
  end;
end;
ys     = a*(x.*v).^alpk.*(l.*u).^alph;
g2     = (z2.*ys).^alpi;
y      = ys.*g2.^alp1;
g1     = z1.*y;
g3     = z3.*y;
xh     = b*(x.*(1-v)).^thetk.*(l.*(1-u)).^theth.*g3.^thet1;
grate  = delh1+xh;
xk     = y-g1-g2-g3-c;
kp     = (delk1*x+xk)./grate;

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

if plotdf;
  plot(x,c),grid, title('consumption'),pause
  plot(x,xk),grid, title('investment in physical capital'),pause
  plot(x,xh),grid, title('investment in human capital'),pause
  plot(x,l),grid, title('hours of work'),pause
  plot(x,grate),grid, title('growth rate'),pause
  plot(x,v),grid, title('fraction of capital in sector 1'),pause
  plot(x,u),grid, title('fraction of labor in sector 1'),pause
  plot(x,kp,x,x),grid, title('capital next period')
end;

%
% Compute steady states.
%

compss = input('Compute steady state? (1=yes, 0=no) ');
if compss;
  policies = [tauxk,tauxh,tauk1,tauk2,tauh1,tauh2,tauc, ...
              z1(1,:)',z2(1,:)',z3(1,:)'];
  istate   = input('Which state for policies? ');
  tem      = find(kp(:,istate)<x(:,1));
  if ~isempty(tem);
    x0     = tem(1);
  else
    x0     = round(ns/2);
  end;
  param    = [a;alpk;alph;b;beta;delk;delh;eta;psi;sigma;thetk;theth;zeta;
              policies(istate,:)'];
  z0       = [x(x0,istate); l(x0,istate); grate(x0,istate); v(x0,istate)];
  z        = secant('twofoc',z0,param);
  kss      = z(1);
  lss      = z(2);
  gamss    = z(3);
  vss      = z(4);
  phi      = tauk11(1,istate)*tauh21(1,istate)*alpk*theth/ ...
             (tauk21(1,istate)*tauh11(1,istate)*alph*thetk);
  uss      = 1/(1-phi+phi/vss);
  yss      = a*(kss*vss)^alpk*(lss*uss)^alph*z2(1,istate)^(1-alpk-alph);
  xhss     = gamss-1+delh;
  xkss     = (gamss-1+delk)*kss;
  css      = (1-z1(1,istate)-z2(1,istate)-z3(1,istate))*yss-xkss;

  disp(' ')
  disp(' ')
  disp(' ')
  disp('Steady state values for Two-sector Growth Model')
  disp('-----------------------------------------------')
  disp(sprintf('   (h''/h-1)*100             =   %7.4f ',(gamss-1)*100))
  disp(sprintf('   c/y                      =   %7.4f ',css/yss))
  disp(sprintf('   g/y                      =   %7.4f ', ...
                                                z1(1,istate)+z2(1,istate)))
  disp(sprintf('   l                        =   %7.4f ',lss))
  disp(sprintf('   xk/y                     =   %7.4f ',xkss/yss))
  disp(sprintf('   xh/y                     =   %7.4f ',xhss/yss))
  disp(sprintf('   k/y                      =   %7.4f ',kss/yss))
  disp(sprintf('   k/h                      =   %7.4f ',kss))
  disp(sprintf('   v                        =   %7.4f ',vss))
  disp(sprintf('   u                        =   %7.4f ',uss))
  disp('-----------------------------------------------')

end;

%
% Simulate time series.
%

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

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

  Ht     = zeros(nc,nt);
  Ct     = zeros(nc,nt);
  Xkt    = zeros(nc,nt);
  Xht    = zeros(nc,nt);
  Gt     = zeros(nc,nt);
  Yt     = zeros(nc,nt);
  Grt    = zeros(nc,nt);
  Kt     = zeros(nc,nt);
  Lt     = zeros(nc,nt);
  St     = zeros(nc,nt);
  Vt     = zeros(nc,nt);
  Ut     = zeros(nc,nt);
  for i=1:nc;
    disp(sprintf('Computing statistics for country %g',i))
    ht(1)  = h0;
    xt(1)  = x0;
    st     = markov(full(pi),nt,i0,i);
    tauxk1t= 1+tauxk(st);
    tauxh1t= 1+tauxh(st);
    tauk11t= 1-tauk1(st);
    tauk21t= 1-tauk2(st);
    tauh11t= 1-tauh1(st);
    tauh21t= 1-tauh2(st);
    tauc1t = 1+tauc(st);
    z1t    = z1(st);
    z2t    = z2(st);
    z3t    = z3(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;
      vt(t,1)    = v(i1,st(t))*bs1 +v(i1+1,st(t))*bs2;
      phi        = tauk11t(t)*tauh21t(t)*alpk*theth/ ...
                  (tauk21t(t)*tauh11t(t)*alph*thetk);
      ut(t,1)    = 1/(1-phi+phi/vt(t));

      tem1       = psi*ct(t);
      tem2       = tauh11t(t)*alph*aa*(xt(t)*vt(t)/ut(t))^alphak* ...
                   z2t(t)^alphai/tauc1t(t);
      lab        = lab0;
      lei        = 1-lab;
      for k=1:10;
        tem      = tem1 - tem2*lab^(-alphak)*lei;
        dtem     = tem2*lab^(-alphak)*(alphak*lei/lab+1);
        lab      = lab - tem/dtem;
        lei      = 1-lab;
      end
      lt(t,1)    = lab;
      ys         = a*(xt(t)*vt(t))^alpk*(lab*ut(t))^alph;
      g2t(t,1)   = (z2t(t)*ys)^alpi;
      yt(t,1)    = ys*g2t(t)^alp1;
      g1t(t,1)   = z1t(t)*yt(t);
      g3t(t,1)   = z3t(t)*yt(t);
      xht(t,1)   = b*(xt(t)*(1-vt(t)))^thetk*(lab*(1-ut(t)))^theth*g3t(t)^thet1;
      gratet(t,1)= delh1+xht(t);
      xkt(t,1)   = yt(t)-g1t(t)-g2t(t)-g3t(t)-ct(t);
      qt(t,1)    = tauk11t(t)*alpk*yt(t)*(1-vt(t))/ ...
                   (tauk21t(t)*thetk*xht(t)*vt(t));
      xt(t+1,1)  = (delk1*xt(t)+xkt(t))/gratet(t);
      ht(t+1,1)  = gratet(t)*ht(t);
    end;
    ht = ht(1:nt);
    xt = xt(1:nt);
    Ht(i,:)  = ht';
    Ct(i,:)  = ct';
    Xkt(i,:) = xkt';
    Xht(i,:) = xht';
    Gt(i,:)  = (g1t+g2t+g3t)';
    Yt(i,:)  = yt';
    Grt(i,:) = gratet';
    Kt(i,:)  = xt';
    Lt(i,:)  = lt';
    St(i,:)  = st';
  end;
  Nt       = (ones(nc,1)*(1+eta/100).^[1:nt]);
  yl       = (Yt+Xht).*Ht./(Lt.*Nt);
  yn       = (Yt+Xht).*Ht./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;


