%SGM1DQ     plots results generated by SGM1DQ.F.

%           Ellen McGrattan, 9-9-94

%
%  Load data written by sgm1dq.f.
%

clear
load sgm1dq.dat
X     = sgm1dq;
clear sgm1dq

%
%  Assign parameters and the consumption function. 
%

alpha = X(1);
beta  = X(2);
delta = X(3);
tau   = X(4);
gamma = X(5);
ptr   = 6;
nx    = X(ptr);
ptr   = ptr+1;
x     = X(ptr:ptr+nx-1);
ptr   = ptr+nx;
ns    = X(ptr);
ptr   = ptr+1;
theta = X(ptr:ptr+ns-1);
ptr   = ptr+ns;
pi    = zeros(ns);
pi(:) = X(ptr:ptr+ns*ns-1);
ptr   = ptr+ns*ns;
n     = nx*ns;
c     = zeros(nx,ns);
c(:)  = X(ptr:ptr+n-1);

%
%  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)]';
else
  xx = input('Mesh for capital (k)      =  ');
end;
nxx = length(xx);

for j=1:ns;
  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;
    bs1    = 0.5*xi*(xi-1);
    bs2    = 1-xi*xi;
    bs3    = 0.5*xi*(xi+1);
    cc(i,j)= c(ii,j)*bs1+c(ii+1,j)*bs2+c(ii+2,j)*bs3;
  end;
end;


%
%  Make pretty plots.
%

s=[];
for i=1:ns;
  str  = num2str(theta(i));
  siz  = length(str);
  s    = [s;[' ',str(1:min(8,siz)),blanks(8-siz)]];
end;
plot(xx,cc),title('Consumption function'),legend(s,-1),pause(2)
plot(xx,xx.^alpha*theta'-cc),title('Investment function'),legend(s,-1)



%
% If test case, report error estimates.
%

if abs(1-delta)< .0001 & abs(1-tau) < .0001;
  disp(' ')
  disp('  ERROR ESTIMATES FOR TEST CASE')
  disp(' ')
  disp('        Theta            Sup-Norm          L2 norm ')
  for j=1:ns
    sum1    = 0;
    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;
      bs1    = 0.5*(xi.*(xi-1));
      bs2    = (1-xi.*xi);
      bs3    = 0.5*(xi.*(xi+1));
      ccc    = c(i,j)*bs1+c(i+1,j)*bs2+c(i+2,j)*bs3;
      d      = (1-beta*alpha)*ax.^alpha*theta(j)-ccc;
      sum1   = sum1+sum(d.*d.*wx);
    end;
    err1(j,1)= sqrt(sum1);
    err2(j,1)= max((1-beta*alpha)*xx.^alpha*theta(j)-cc(:,j));
    disp(sprintf('   %12.6f      %12.6f      %12.6f ',theta(j),err2(j),err1(j)))
  end;
end;



