%SGM1DL    plot results from SGM1DL.F.
%

%          Ellen McGrattan, 9-7-94


%
% Load data written by sgm1dl.f.
%

load sgm1dl.dat
X     = sgm1dl;
clear sgm1dl

%
% 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 = exp(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);
p     = pi';
Pi    = (eye(ns-1)-p(1:ns-1,1:ns-1)+p(1:ns-1,ns)*ones(1,ns-1))\p(1:ns-1,ns);
Pi    = [Pi;1-sum(Pi)];

%
% Construct decision functions for investment and capital.
%

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

%
% Compute the steady state for capital.
%

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

%
% 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;
legend off
if ns>1; mesh(x*ones(1,ns),ones(nx,1)*theta',c), ...
   title('Consumption function'),pause(2);end;  
plot(x,c,x,c,'x'),title('Consumption function'),grid,legend(s,-1), pause(2)
plot(x,in,x,in,'x'),title('Investment function'),grid,legend(s,-1),pause(2)
str  = [' steady  ';
        ' state   '];
s    = [s;str];
plot(x,kp,sk,sk,'o'),title('Capital function'),grid,legend(s,-1) 

%
% If test case, check errors.
%

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;
    best = 0;
    for i=1:nx-1;
      x1    = x(i);
      x2    = x(i+1);
      c1    = c(i,j);
      c2    = c(i+1,j);
      dx    = x2-x1;
      [a,w] = qgausl(x1,x2,20);
      xp    = (x2-a)/dx;
      xm    = (a-x1)/dx;
      d     = (1-beta*alpha)*a.^alpha*theta(j)-c1*xp-c2*xm;
      sum1  = sum1+sum(d.*d.*w);
      xmax  = ((c2-c1)/(x2-x1)/(alpha*theta(j)*(1-alpha*beta)))^(1/(alpha-1));
      dmax  = (1-alpha*beta)*theta(j)*xmax^alpha-(c1*(x2-xmax)+ ...
              c2*(xmax-x1))/dx;
      if xmax>x2 | xmax <x1; dmax=0; end;
      best  = max([best;dmax;abs((1-alpha*beta)*theta(j)*x1^alpha-c1);
                   abs((1-alpha*beta)*theta(j)*x2^alpha-c2)]);
    end;
    err2(j,1) = sqrt(sum1);
    err1(j,1) = best;
    disp(sprintf('   %12.6f      %12.6f      %12.6f ',theta(j),err1(j),err2(j)))
  end; 
end;

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

nt    = 300;
i0    = find(rand<cumsum(Pi));
i0    = i0(1);
st    = markov(pi,nt,i0);
xt(1) = ((1-beta*(1-delta))/(theta(st(1))*beta*alpha))^(1/(alpha-1));

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;
  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;
  Y          = theta(st(t))*xt(t)^alpha;
  xt(t+1,1)  = Y+(1-delta)*xt(t)-ct(t);
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.*theta(st)),xlabel('time'),title('Output'),grid

