%CU3D    is a script file that analyzes results of CU3DA.f.
%        cu3da.f computes decision functions for the stochastic 
%        growth model with government capital and varying capacity
%        utilization: 
%
%              max           E[ sum  beta^t (1-a[t])*{n[t]*u(c1[t],l[t]) 
%         c1[t],i[t],h[t],       t             +(1-n[t]*u(c0[t],1) }|x[0]]
%           c0[t],n[t]
%
%              subject to
%
%                  c[t]+i[t] = lam[t]*f(k[t],kg[t],n[t],h[t])-ig[t]-cg[t]
%                  c[t]      = n[t]*c1[t]+(1-n[t])*c0[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,cg,army:  discrete Markov processes
%

%        Ellen McGrattan, 4-16-93
%        Revised, ERM, 7-19-93

clear

% read in answer:

load cu3d.dat
X = cu3d;
clear cu3d

beta  = X(1);
delta = X(2);
gamma = X(3);
lam   = X(4);
omega = X(5);
pa    = X(6);
pb    = X(7);
phi   = X(8);
psi   = X(9);
rho   = X(10);
sig   = X(11);
theta = X(12);
zeta  = X(13);
ptr   = 14;
ns    = X(ptr);
ptr   = ptr+1;
ig    = X(ptr:3:ptr-1+3*ns);
cg    = 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;
c1    = zeros(n,ns);
c1(:) = X(ptr:ptr-1+n*ns);
delt1 = 1.d0-delta;
gamm1 = 1.d0-gamma;
ome1  = 1.d0-omega;
rho1  = 1.d0-rho;
thet1 = 1.d0-theta;
cwgt  = gamma*ome1;
cwgt1 = cwgt-1.d0;
hwgt  = gamm1*ome1;
hwgt1 = hwgt-1.d0;
Pi    = ones(ns,1)/ns;
for i = 1:20;
  Pi  = (Pi'*pi)';
end;



% for each state, compute the hours (h), the fraction working (n0), 
% and consumption for those not working (c0)

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)*cg';
aa    = ones(n,1)*army';

h     = 0.5*ones(n,ns);
n0    = 0.1*ones(n,ns);
zkk   = lam*zz.*(phi*xx.^rho+(1-phi)*yy.^rho).^(theta/rho);
     
for  m=1:30
  c0   = c1.*(1-h).^(hwgt/cwgt1);
% dpn  = pa+2*pb*n0;
% d2pn = 2*pb;
% dpn  = 3*pb*n0.*n0;
% d2pn = 6*pb*n0;
  tem  = (n0-pa)*pb;
  dpn  = 3*tem.*tem.*(1+sign(tem))*pb;
  d2pn = 6*tem.*(1+sign(tem))*pb*pb;
  fj11 =-d2pn./(c1.^cwgt)-gamma*thet1*theta*(1-h).^hwgt./c1.*zkk.*n0.^ ...
         (-theta-1).*h;
  fj22 = zkk;
  fj12 = hwgt/ome1.*(1-h).^hwgt1-gamma*hwgt*(1-h).^hwgt1./c1.*(c1-c0- ...
         thet1*zkk.*n0.^(-theta).*h)-gamma*thet1*(1-h).^hwgt./c1.*zkk.* ...
         n0.^(-theta);
  fj21 = -gamm1*theta*n0.^(-thet1).*c1/gamma;
  det  = fj11.*fj22-fj12.*fj21;
  f1   = gamm1/gamma*n0.^theta.*c1-zkk.*(1-h);
  f2   = (1-h).^hwgt/ome1-(c0./c1).^cwgt/ome1-dpn./(c1.^cwgt)-gamma* ...
         (1-h).^hwgt./c1.*(c1-c0-thet1*zkk.*n0.^(-theta).*h);
  h    =  h - (fj11.*f1+fj21.*f2)./det;
  n0   =  n0 - (fj12.*f1+fj22.*f2)./det;
end;


% steady states

zss  = 1;
igss = ig(1);
cgss = cg(1);
ass  = army(1);
kgss = igss/delta;
par  = [X(1:12);igss;cgss;ass];
ss   = secant('sscu',[.6;.4;.2;.8],par);
hss  = ss(1);
nss  = ss(2);
c1ss = ss(3);
kss  = ss(4);

% simulations:

T       = 1000;
xt(1)   = kss;
yt(1)   = kgss;
zt(1)   = zss;
igt(1)  = igss;
cgt(1)  = cgss;
at(1)   = ass;

randn('seed',0);
et      = sig*randn(T,1);
st      = markov(pi,T);
st      = st(:);
nsim    = 1;
%load simil
%[T,nsim]= size(simil);
%st     = simil;
xt      = xt*ones(1,nsim);
yt      = yt*ones(1,nsim);
zt      = zt*ones(1,nsim);
igt     = igt*ones(1,nsim);
cgt     = cgt*ones(1,nsim);
at      = at*ones(1,nsim);

for j=1:nsim;
  for t=1:T;
    l        = st(t,j);
    igt(t,j) = ig(l);
    cgt(t,j) = cg(l)+zeta*log(zt(t,j));
    at(t,j)  = army(l);

    low      = find(xt(t,j)>x);
    high     = find(xt(t,j)<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,j)>y);
    high     = find(yt(t,j)<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);
    jj       = low;
    low      = find(zt(t,j)>z);
    high     = find(zt(t,j)<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+(jj-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,j))/(x2-x1);
    xm       = (xt(t,j)-x1)/(x2-x1);
    yp       = (y2-yt(t,j))/(y2-y1);
    ym       = (yt(t,j)-y1)/(y2-y1);
    zp       = (z2-zt(t,j))/(z2-z1);
    zm       = (zt(t,j)-z1)/(z2-z1);
    c1t(t,j) = c1(nd1,l)*xp*yp*zp+c1(nd2,l)*xm*yp*zp+c1(nd3,l)*xm*ym*zp+ ...
               c1(nd4,l)*xp*ym*zp+c1(nd5,l)*xp*yp*zm+c1(nd6,l)*xm*yp*zm+ ...
               c1(nd7,l)*xm*ym*zm+c1(nd8,l)*xp*ym*zm;
    ht(t,j)  = 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;
    nt(t,j)  = n0(nd1,l)*xp*yp*zp+n0(nd2,l)*xm*yp*zp+n0(nd3,l)*xm*ym*zp+ ...
               n0(nd4,l)*xp*ym*zp+n0(nd5,l)*xp*yp*zm+n0(nd6,l)*xm*yp*zm+ ...
               n0(nd7,l)*xm*ym*zm+n0(nd8,l)*xp*ym*zm;
    c0t(t,j) = c0(nd1,l)*xp*yp*zp+c0(nd2,l)*xm*yp*zp+c0(nd3,l)*xm*ym*zp+ ...
               c0(nd4,l)*xp*ym*zp+c0(nd5,l)*xp*yp*zm+c0(nd6,l)*xm*yp*zm+ ...
               c0(nd7,l)*xm*ym*zm+c0(nd8,l)*xp*ym*zm;
    ft(t,j)  = zt(t,j)*lam*(phi*xt(t,j)^rho+(1-phi)*yt(t,j)^rho)^ ...
               (theta/rho)*nt(t,j)^(1-theta)*ht(t,j);
    xt(t+1,j)= ft(t,j)+(1-delta)*xt(t,j)-nt(t,j)*c1t(t,j)-(1-nt(t,j))* ...
	       c0t(t,j)-igt(t,j)-cgt(t,j);
    yt(t+1,j)= (1-delta)*yt(t,j)+igt(t,j);
    zt(t+1,j)= exp(psi*log(zt(t,j))+et(t));
  end;

  it(:,j)  = xt(2:T+1,j)-(1-delta)*xt(1:T,j);
  xt       = xt(1:T,:);
  yt       = yt(1:T,:);
  zt       = zt(1:T,:);
  Ht(:,j)  = nt(:,j).*ht(:,j);
  prt(:,j) = ft(:,j)./Ht(:,j);
  gt(:,j)  = cgt(:,j)+igt(:,j);
  ct(:,j)  = nt(:,j).*c1t(:,j)+(1-nt(:,j)).*c0t(:,j);

  t   = ptime(1900,T,1);
  tem = [ft(:,j),ct(:,j),it(:,j),gt(:,j),ht(:,j),nt(:,j),prt(:,j)];
  tem = tem(47:T,:);
  stat(j,:) = 100*std(log(tem)-hptrend(log(tem),1600));
  disp(' ')
  disp(' ')
  disp(sprintf('std. dev. of output               = %g ',stat(j,1)))
  disp(sprintf('std. dev. of consumption          = %g ',stat(j,2)))
  disp(sprintf('std. dev. of investment           = %g ',stat(j,3)))
  disp(sprintf('std. dev. of government spending  = %g ',stat(j,4)))
  disp(sprintf('std. dev. of hours/worker         = %g ',stat(j,5)))
  disp(sprintf('std. dev. of fraction employed    = %g ',stat(j,6)))
  disp(sprintf('std. dev. of productivity         = %g ',stat(j,7)))
  disp(' ')
  disp(' ')
  tem = corrcoef([ft(2:T-1,j),ht(1:T-2,j),ht(2:T-1,j),ht(3:T,j)]);
  disp(sprintf('correlation of y[t] and h[t-1],h[t],h[t-1]: %g,%g,%g', ...
                tem(1,2),tem(1,3),tem(1,4)))
  tem = corrcoef([ft(2:T-1,j),nt(1:T-2,j),nt(2:T-1,j),nt(3:T,j)]);
  disp(sprintf('correlation of y[t] and n[t-1],n[t],n[t-1]: %g,%g,%g', ...
                tem(1,2),tem(1,3),tem(1,4)))
  disp(' ')
  disp(' ')
  disp(' ')

  plot(t,ht(:,j),t,nt(:,j)),pause

end;

return
plot(t,ft,t,ct,t,it,t,gt,t,igt),grid
pause
plot(t,ct./ft,t,it./ft,t,cgt./ft,t,igt./ft),grid
pause
plot(t,xt./ft,t,yt./ft),grid
pause
plot(t,gt./ft,t,cgt./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



