%SGM2DL   plots results from SGM2DL.F
%         

%         Ellen McGrattan, 9-9-94

%
%  Load data written by sgm2dl.f.
%

clear

load sgm2dl.dat
X     = sgm2dl;
clear sgm2dl

%
%  Assign parameters and the consumption function. 
%

alpha = X(1);
beta  = X(2);
delta = X(3);
rho   = X(4);
sige  = X(5);
tau   = X(6);
ptr   = 7;
nx    = X(ptr);
ptr   = ptr+1;
x     = X(ptr:ptr+nx-1);
ptr   = ptr+nx;
ny    = X(ptr);
ptr   = ptr+1;
y     = X(ptr:ptr+ny-1);
ptr   = ptr+ny;
n     = nx*ny;
c     = zeros(nx,ny);
c(:)  = X(ptr:ptr+n-1);
theta = exp(0.5*log((1+y)./(1-y)));

%
% Construct decision functions for investment and capital.
%

in    = x.^alpha*theta'-c;
kp    = x.^alpha*theta'+(1-delta)*x*ones(1,ny)-c;

%
% Compute the steady state for capital.
%

sk    = ((1-beta*(1-delta))/(beta*alpha))^(1/(alpha-1));

%
% Make pretty plots.
%

s=[];
for i=1:ny;
  str  = num2str(theta(i));
  siz  = length(str); 
  s    = [s;[' ',str(1:min(8,siz)),blanks(8-siz)]];
end;
legend off
if ny>1;
  mesh(x*ones(1,ny),ones(nx,1)*theta',c),title('Consumption function'),pause(2)
end;
plot(x,c,x,c,'+'),title('Consumption function'),legend(s,-1),grid,pause(2)
plot(x,in,x,in,'+'),title('Investment function'),legend(s,-1),grid,pause(2)
str  = [' steady  ';
        ' state   '];
s    = [s;str];
plot(x,kp,sk,sk,'o'),title('Capital function'),grid,legend(s,-1),grid,pause(2)

%
% If test case, report error estimates.
%

if abs(1-delta)< .0001 & abs(1-tau) < .0001;
  sum1    = 0;
  for j=1:ny-1;
    y1      = y(j);
    y2      = y(j+1);
    dy      = y2-y1;
    [ay,wy] = qgausl(y1,y2,20);
    yp      = (y2-ay)/dy;
    ym      = (ay-y1)/dy;
    for i=1:nx-1;
      x1      = x(i);
      x2      = x(i+1);
      dx      = x2-x1;
      [ax,wx] = qgausl(x1,x2,20);
      xp      = (x2-ax)/dx;
      xm      = (ax-x1)/dx;
      nd1     = i+(j-1)*nx; 
      nd2     = nd1+1;
      nd4     = nd1+nx;
      nd3     = nd4+1;
      cc      = c(nd1)*xp*yp'+c(nd2)*xm*yp'+c(nd3)*xm*ym'+c(nd4)*xp*ym';
      d       = (1-beta*alpha)*ax.^alpha*sqrt((1+ay)./(1-ay))'-cc;
      sum1    = sum1+sum(sum(d.*d.*(wx*wy')));
    end;
  end;
  err1 = sqrt(sum1);
  disp(' ')
  disp('  ERROR ESTIMATES FOR TEST CASE')
  disp(' ')
  disp(sprintf('  L2  - norm           %12.6f',err1))
  disp('  Sup - norm           ')
  err2  = 0;
  for j=1:ny;
    for i=1:nx;
      d    = (1-beta*alpha)*x(i)^alpha*sqrt((1+y(j))/(1-y(j)))-c(i,j);
      err2 = max(err2,abs(d));
    end;
  end;
  disp(sprintf('    w/ original grid   %12.6f',err2))
  xx = expstret(.00005,x(1),x(nx),100,1);
  yy = [y(1):.01*(y(ny)-y(1)):y(ny)]';
  err3  = 0;
  for j=1:length(yy);
    jj     = find(y<yy(j));
    if (jj==[]);  jj=1;   end;
    jj     = jj(length(jj));
    if (jj>ny-1); jj=ny-1; end;
    y1     = y(jj);
    y2     = y(jj+1);
    dy     = y2-y1;
    yp     = (y2-yy(j))/dy;
    ym     = (yy(j)-y1)/dy;
    for i=1:length(xx);
      ii     = find(x<xx(i));
      if (ii==[]);  ii=1;   end;
      ii     = ii(length(ii));
      if (ii>nx-1); ii=nx-1; end;
      x1     = x(ii);
      x2     = x(ii+1);
      dx     = x2-x1;
      xp     = (x2-xx(i))/dx;
      xm     = (xx(i)-x1)/dx;
      nd1    = ii+(jj-1)*nx;
      nd2    = nd1+1;
      nd4    = nd1+nx;
      nd3    = nd4+1;
      cc(i,j)= c(nd1)*xp*yp+c(nd2)*xm*yp+c(nd3)*xm*ym+c(nd4)*xp*ym;
      d      = (1-beta*alpha)*xx(i)^alpha*exp(.5*log(1+yy(j))- ...
               .5*log(1-yy(j)))'-cc(i,j);
      err3   = max(err3,abs(d));
    end;
  end;
  disp(sprintf('    w/ tight grid      %12.6f',err3))
  disp(' ')
end;

%
% Generate time series for consumption, investment, and output and plot it.
%

nt    = 1000;
xt(1) = sk;
yt(1) = 0;

for t=1:nt;
  i1         = find(x<xt(t));
  if (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;
  j1         = find(y<yt(t));
  if (j1==[]);  j1=1;   end;
  j1         = j1(length(j1));
  if (j1>ny-1); j1=ny-1; end;
  y1         = y(j1);
  y2         = y(j1+1);
  dy         = y2-y1;
  xi         = 2*(xt(t)-x1)/dx-1;
  eta        = 2*(yt(t)-y1)/dy-1;
  bs1        = 0.25*(1-xi)*(1-eta);
  bs2        = 0.25*(1+xi)*(1-eta);
  bs3        = 0.25*(1+xi)*(1+eta);
  bs4        = 0.25*(1-xi)*(1+eta);
  nd1        = (j1-1)*nx+i1;
  nd2        = nd1+1;
  nd4        = nd1+nx;
  nd3        = nd4+1;
  ct(t,1)    = c(nd1)*bs1+c(nd2)*bs2+c(nd3)*bs3+c(nd4)*bs4;
  tech       = 0.5*(log((1+yt(t))/(1-yt(t))));
  Y          = exp(tech)*xt(t)^alpha;
  xt(t+1,1)  = Y+(1-delta)*xt(t)-ct(t);
  yt(t+1,1)  = tanh(rho*tech+sige*randn);
end;
legend off
plot(1:nt,xt(1:nt)),xlabel('time'),title('Capital stock'),grid,pause(1)
plot(1:nt,ct),xlabel('time'),title('Consumption'),grid,pause(1)
plot(1:nt,xt(2:nt+1)-(1-delta)*xt(1:nt)),xlabel('time'), ...
  title('Investment'),grid,pause(1)
plot(1:nt,xt(1:nt).^alpha.*sqrt((1+yt(1:nt))./(1-yt(1:nt)))), ...
   xlabel('time'),title('Output'),grid

