%SGM2DQ     plots results generated by SGM2DQ.F.

%           Ellen McGrattan, 9-9-94

%
%  Load data written by sgm2dq.f.
%

clear
load sgm2dq.dat
X     = sgm2dq;
clear sgm2dq

%
%  Assign parameters and the consumption function. 
%

alpha = X(1);
beta  = X(2);
delta = X(3);
rho   = X(4);
sige  = X(5);
tau   = X(6);
gamma = X(7);
ptr   = 8;
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     = X(ptr:ptr+n-1);
theta = sqrt((1+y)./(1-y));

%
%  Compute the steady state for capital.
%

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


%
% Refine grid for plotting.
%

mtype = input('Do you want an evenly spaced mesh? (1=yes) ');
if mtype==1;
  dif = input('What mesh size for capital? ');
  xx  = [x(1):dif:x(nx),x(nx)]';
  dif = input('What mesh size for technology? ');
  yy  = [y(1):dif:y(ny),y(ny)]';
  tt  = sqrt((1+yy)./(1-yy));
else
  xx = input('Mesh for capital (k)      =  ');
  yy = input('Mesh for tanh(log(theta)) =  ');
  tt = sqrt((1+yy)./(1-yy));
end;
nxx = length(xx);
nyy = length(yy);

for j=1:nyy;
  jj     = find(y<yy(j));
  if (jj==[]);  jj=1;   end;
  jj     = jj(length(jj));
  if (jj>ny-1); jj=ny-1; end;
  jj     = jj-rem(jj+1,2);
  y1     = y(jj);
  y2     = y(jj+2);
  for i=1:nxx;
    ii     = find(x<xx(i));
    if (ii==[]);  ii=1;   end;
    ii     = ii(length(ii));
    if (ii>nx-1); ii=nx-1; end;
    ii     = ii-rem(ii+1,2);
    x1     = x(ii);
    x2     = x(ii+2);
    xi     = 2*(xx(i)-x1)/(x2-x1)-1;
    eta    = 2*(yy(j)-y1)/(y2-y1)-1;
    bs1    = 0.25*xi*(xi-1)*eta*(eta-1);
    bs2    = 0.25*xi*(xi+1)*eta*(eta-1);
    bs3    = 0.25*xi*(xi+1)*eta*(eta+1);
    bs4    = 0.25*xi*(xi-1)*eta*(eta+1);
    bs5    = 0.5*(1-xi*xi)*eta*(eta-1);
    bs6    = 0.5*xi*(xi+1)*(1-eta*eta);
    bs7    = 0.5*(1-xi*xi)*eta*(eta+1);
    bs8    = 0.5*xi*(xi-1)*(1-eta*eta);
    bs9    = (1-xi*xi)*(1-eta*eta);
    nd1    = ii+nx*(jj-1);
    nd5    = nd1+1;
    nd2    = nd1+2;
    nd4    = nd1+2*nx;
    nd7    = nd4+1;
    nd3    = nd4+2;
    nd8    = nd1+nx;
    nd9    = nd8+1;
    nd6    = nd8+2;
    cc(i,j)= c(nd1)*bs1+c(nd2)*bs2+c(nd3)*bs3+c(nd4)*bs4+ ...
             c(nd5)*bs5+c(nd6)*bs6+c(nd7)*bs7+c(nd8)*bs8+ ...
             c(nd9)*bs9;
  end;
end;


%
%  Make pretty plots.
%

legend off
s=[];
for i=1:nyy;
  str  = num2str(tt(i));
  siz  = length(str);
  s    = [s;[' ',str(1:min(8,siz)),blanks(8-siz)]];
end;
if nyy>1;
  mesh(xx*ones(1,nyy),ones(nxx,1)*yy',cc),title('Consumption function'),pause(2)
end;
plot(xx,cc),title('Consumption function'), pause(2)
plot(xx,xx.^alpha*tt'-cc),title('Investment function')


%
% If test case, report error estimates.
%


if abs(1-delta)< .0001 & abs(1-tau) < .0001;
  disp(' ')
  disp('  ERROR ESTIMATES FOR TEST CASE')
  disp(' ')
  sum1    = 0;
  for j=1:2:ny-2;
    y1      = y(j);
    y2      = y(j+2);
    [ay,wy] = qgausl(y1,y2,20);
    for i=1:2:nx-2;
      x1    = x(i);
      x2    = x(i+2);
      [ax,wx] = qgausl(x1,x2,20);
      xi     = 2*(ax-x1)/(x2-x1)-1;
      eta    = 2*(ay-y1)/(y2-y1)-1;
      bs1    = 0.25*(xi.*(xi-1))*(eta.*(eta-1))';
      bs2    = 0.25*(xi.*(xi+1))*(eta.*(eta-1))';
      bs3    = 0.25*(xi.*(xi+1))*(eta.*(eta+1))';
      bs4    = 0.25*(xi.*(xi-1))*(eta.*(eta+1))';
      bs5    = 0.5*(1-xi.*xi)*(eta.*(eta-1))';
      bs6    = 0.5*(xi.*(xi+1))*((1-eta.*eta))';
      bs7    = 0.5*(1-xi.*xi)*(eta.*(eta+1))';
      bs8    = 0.5*(xi.*(xi-1))*(1-eta.*eta)';
      bs9    = (1-xi.*xi)*(1-eta.*eta)';
      nd1    = i+nx*(j-1);
      nd5    = nd1+1;
      nd2    = nd1+2;
      nd4    = nd1+2*nx;
      nd7    = nd4+1;
      nd3    = nd4+2;
      nd8    = nd1+nx;
      nd9    = nd8+1;
      nd6    = nd8+2;
      ccc    = c(nd1)*bs1+c(nd2)*bs2+c(nd3)*bs3+c(nd4)*bs4+ ...
               c(nd5)*bs5+c(nd6)*bs6+c(nd7)*bs7+c(nd8)*bs8+ ...
               c(nd9)*bs9;
      d       = (1-beta*alpha)*ax.^alpha*(sqrt((1+ay)./(1-ay)))'-ccc;
      sum1    = sum1+sum(sum(d.*d.*(wx*wy')));
    end;
  end;
  err1 = sqrt(sum1);
  disp(sprintf('  L2  - norm           %12.6f',err1))
  err2 = max(max((1-beta*alpha)*xx.^alpha*(sqrt((1+yy)./(1-yy)))'-cc));
  disp(sprintf('  Sup - norm           %12.6f',err2))
end;

