%WAR3D   is a script file that analyzes results of WAR3DA.f/WAR3DB.f.
%        War3d*.f computes decision functions for the stochastic 
%        growth model with government capital: 
%
%            max           E[ sum  beta^t u(c[t],l[t]) | k[0],kg[0],lam[0]]
%         c[t],i[t],h[t]       t
%
%              subject to
%
%                        c[t]+i[t] = lam[t]*f(k[t],kg[t],h[t])-ig[t]-bg[t]
%                        k[t+1]    = (1-delta)*k[t]+i[t]
%                        kg[t+1]   = (1-delta)*kg[t]+ig[t]
%                        lam[t+1]  = lam[t]^psi*eps[t]
%                        ig,bg,army:  discrete Markov processes
%

%        Ellen McGrattan, 3-5-93


% read in answer:

load war3d.dat
X = war3d;
clear war3d

alpha = X(1);
beta  = X(2);
delta = X(3);
lam   = X(4);
phi   = X(5);
psi   = X(6);
rho   = X(7);
sig   = X(8);
theta = X(9);
ptr   = 10;
ns    = X(ptr);
ptr   = ptr+1;
ig    = X(ptr:3:ptr-1+3*ns);
bg    = X(ptr+1:3:ptr-1+3*ns);
army  = X(ptr+2:3:ptr-1+3*ns);
ptr   = ptr+3*ns;
pi    = zeros(ns);
pi(:) = X(ptr:ptr-1+ns*ns);
ptr   = ptr+ns*ns;
nx1   = X(ptr); 
ptr   = ptr+1;
ny1   = X(ptr); 
ptr   = ptr+1;
nz1   = X(ptr); 
ptr   = ptr+1;
x     = X(ptr:ptr-1+nx1);
ptr   = ptr+nx1;
y     = X(ptr:ptr-1+ny1);
ptr   = ptr+ny1;
z     = X(ptr:ptr-1+nz1);
ptr   = ptr+nz1;
n     = nx1*ny1*nz1;
c     = zeros(n,ns);
c(:)  = X(ptr:ptr-1+n*ns);
thet1 = (1-theta);
Pi    = ones(ns,1)/ns;
for i = 1:20;
  Pi = (Pi'*pi)';
end;


% for each state, compute the hours decision

xx    = kron(ones(ny1*nz1,1),x)*ones(1,ns);
yy    = kron(ones(nz1,1),kron(y,ones(nx1,1)))*ones(1,ns);
zz    = kron(z,ones(nx1*ny1,1))*ones(1,ns);
ii    = ones(n,1)*ig';
bb    = ones(n,1)*bg';
aa    = ones(n,1)*army';
cnst  = thet1*lam*zz.*(phi*xx.^rho+(1-phi)*yy.^rho).^(theta/rho)./c;
h     = .1*ones(c);
for i=1:5;
  h   = (cnst.*(1-aa)-alpha*thet1*h.^theta)./(cnst+alpha*theta*h.^(-thet1));
end;

% steady states

zss   = 1;
igss  = ig(1);
bgss  = bg(1);
ass   = army(1);
kgss  = igss/delta;
y0    = [.1;.3];
del   = .0001*eye(2);
for j=1:10;
 for i=1:2;
    yp = y0 + del(:,i);
    k0 = y0(1);
    h0 = y0(2);
    kp = yp(1);
    hp = yp(2);
    c0 = zss*lam*(phi*k0^rho+(1-phi)*kgss^rho)^(theta/rho)*h0^(1-theta)- ...
         delta*k0-igss-bgss;
    cp = zss*lam*(phi*kp^rho+(1-phi)*kgss^rho)^(theta/rho)*hp^(1-theta)- ...
         delta*kp-igss-bgss;
    w0 = zss*lam*thet1*(phi*k0^rho+(1-phi)*kgss^rho)^(theta/rho)*h0^(-theta);
    wp = zss*lam*thet1*(phi*kp^rho+(1-phi)*kgss^rho)^(theta/rho)*hp^(-theta);
    r0 = zss*lam*theta*phi*k0^(rho-1)*(phi*k0^rho+(1-phi)*kgss^rho)^ ...
         (theta/rho-1)*h0^(1-theta);
    rp = zss*lam*theta*phi*kp^(rho-1)*(phi*kp^rho+(1-phi)*kgss^rho)^ ...
         (theta/rho-1)*hp^(1-theta);
    f0 = [alpha*c0-w0*(1-h0-ass);
          beta*r0-1+beta*(1-delta)];
    fp = [alpha*cp-wp*(1-hp-ass);
          beta*rp-1+beta*(1-delta)];
    fj(:,i) = (fp-f0)/del(i,i);
  end;
  y1 = y0-fj\f0;
  y0 = y1;
end;
kss = y1(1);
hss = y1(2);
css = c0;

% simulations:

T     = 100;
xt(1) = kss;
yt(1) = kgss;
zt(1) = zss;
igt(1)= igss;
bgt(1)= bgss;
at(1) = ass;
rand('normal');
rand('seed',3497);
et    = sig*rand(T,1);
st    = markov(pi,T);

for t=1:T;
  l        = st(t);
  igt(t,1) = ig(l);
  bgt(t,1) = bg(l);
  at(t,1)  = army(l);

  low      = find(xt(t)>x);
  high     = find(xt(t)<x);
  if isempty(low); 
    low  = 1; 
    high = 2; 
  elseif isempty(high)
    low  = nx1-1; 
    high = nx1; 
  else
    low  = low(length(low));
    high = low+1;
  end;
  x1       = x(low);
  x2       = x(high);
  i        = low;
  low      = find(yt(t)>y);
  high     = find(yt(t)<y);
  if isempty(low); 
    low  = 1; 
    high = 2; 
  elseif isempty(high)
    low  = ny1-1; 
    high = ny1; 
  else
    low  = low(length(low));
    high = low+1;
  end;
  y1       = y(low);
  y2       = y(high);
  j        = low;
  low      = find(zt(t)>z);
  high     = find(zt(t)<z);
  if isempty(low); 
    low  = 1; 
    high = 2; 
  elseif isempty(high)
    low  = nz1-1; 
    high = nz1; 
  else
    low  = low(length(low));
    high = low+1;
  end;
  if isempty(high); 
    high = nz1; 
  else 
    high = high(1);
  end;
  z1       = z(low);
  z2       = z(high);
  k        = low;
  nd1      = (k-1)*nx1*ny1+(j-1)*nx1+i;
  nd2      = nd1+1;
  nd4      = nd1+nx1;
  nd3      = nd4+1;
  nd5      = nd1+nx1*ny1;
  nd6      = nd5+1;
  nd8      = nd5+nx1;
  nd7      = nd8+1;
  xp       = (x2-xt(t))/(x2-x1);
  xm       = (xt(t)-x1)/(x2-x1);
  yp       = (y2-yt(t))/(y2-y1);
  ym       = (yt(t)-y1)/(y2-y1);
  zp       = (z2-zt(t))/(z2-z1);
  zm       = (zt(t)-z1)/(z2-z1);
  ct(t,1)  = c(nd1,l)*xp*yp*zp+c(nd2,l)*xm*yp*zp+c(nd3,l)*xm*ym*zp+ ...
             c(nd4,l)*xp*ym*zp+c(nd5,l)*xp*yp*zm+c(nd6,l)*xm*yp*zm+ ...
             c(nd7,l)*xm*ym*zm+c(nd8,l)*xp*ym*zm;
  ht(t,1)  = h(nd1,l)*xp*yp*zp+h(nd2,l)*xm*yp*zp+h(nd3,l)*xm*ym*zp+ ...
             h(nd4,l)*xp*ym*zp+h(nd5,l)*xp*yp*zm+h(nd6,l)*xm*yp*zm+ ...
             h(nd7,l)*xm*ym*zm+h(nd8,l)*xp*ym*zm;
  %or ...
  H       = .5;
  cnst    = thet1*zt(t)*lam*(phi*xt(t)^rho+(1-phi)*yt(t)^rho)^(theta/rho)/ct(t);
  H       = (cnst*(1-at(t))-alpha*thet1*H^theta)/(cnst+alpha*theta*H^(-thet1));
  H       = (cnst*(1-at(t))-alpha*thet1*H^theta)/(cnst+alpha*theta*H^(-thet1));
  H       = (cnst*(1-at(t))-alpha*thet1*H^theta)/(cnst+alpha*theta*H^(-thet1));
  H       = (cnst*(1-at(t))-alpha*thet1*H^theta)/(cnst+alpha*theta*H^(-thet1));
  H       = (cnst*(1-at(t))-alpha*thet1*H^theta)/(cnst+alpha*theta*H^(-thet1));
  ht(t)   = H;
  ft(t,1) = zt(t)*lam*(phi*xt(t)^rho+(1-phi)*yt(t)^rho)^(theta/rho)* ..
              ht(t)^(1-theta);
  xt(t+1,1)= ft(t)+(1-delta)*xt(t)-ct(t)-igt(t)-bgt(t);
  yt(t+1,1)= (1-delta)*yt(t)+igt(t);
  zt(t+1,1)= exp(psi*log(zt(t))+et(t));
end;

prt= ft./ht;
it = xt(2:T+1)-(1-delta)*xt(1:T);
xt = xt(1:T);
yt = yt(1:T);
zt = zt(1:T);
gt = bgt+igt;

t = [1910:1909+T]';
plot(t,ft,t,ct,t,it,t,gt,t,igt),grid
pause
plot(t,ct./ft,t,it./ft,t,bgt./ft,t,igt./ft),grid
pause
plot(t,xt./ft,t,yt./ft),grid
pause
plot(t,gt./ft,t,bgt./ft,t,igt./ft,'-.'),grid
pause
plot(t,prt/mean(prt),t,ft/mean(ft),t,ht/mean(ht),'-.'),grid
pause
plot(t,at./ht),grid

