%TABLEA    This file is used to compute the results of Table 1, column 2,
%          Table 2, column 2, and Table 3, column 2 in McGrattan, ``A 
%          Progress Report on Business Cycle Models.''
%
%          The optimization problem being solved is: 
%
%              v(k,z) =  max  {u(c+pi*g,l)+beta E[v(k',z')|k,z]}
%                      c,i,l,k'
%
%                s.t.  c+i    = r*k+w*n -tau*(r-delta)*k-phi*w*n+Tr
%                      k'     = (1-delta)*k + i
%                      z'     = rho0 + rho*z + eps,  std(eps) = sig
%                      u(c,l) = ln(c) + gamma*ln(l)
%                      f(k,n) = k^theta*n^(1-theta)
%
%                      solution:   i(k,z),n(k,z)
 
%          Ellen McGrattan, 12-21-92
%          Revised, 11-22-94, ERM

clear 

% Parameters

case = input('Which table? (1, 2, or 3)  ');
disp(' ')

sig        = zeros(4);
if case==1;
  theta    = 0.355;
  delta    = 0.023;
  beta     = 0.99;
  gamma    = 2.36;
  pi       = 0;
  rho      = [0.9739,0,0,0;
              0,1e-10,0,0;
              0,0,1e-10,0;
              0,0,0,1e-10];
  rho0     = [1.0*(1-rho(1,1)),0*(1-rho(2,2)),0*(1-rho(3,3)),0*(1-rho(4,4))]';
  sig      = zeros(4);
  sig(1,1) = 0.0096;
elseif case==2;
  theta    = 0.359;
  delta    = 0.023;
  beta     = 0.99;
  gamma    = 2.33;
  pi       = 0;
  rho      = [0.9739,0,0,0;
              0,.9687,0,0;
              0,0,.9755,0;
              0,0,0,.9703];
  rho0     = [1.0*(1-rho(1,1)),.195*(1-rho(2,2)),.5*(1-rho(3,3)), ....
              .23*(1-rho(4,4))]';
  sig(1,1) = 0.0096;
  sig(2,1) = 0.0033;
  sig(2,2) = 0.0061;
  sig(3,1) = 0.0012;
  sig(3,2) = 0.0003;
  sig(3,3) = 0.0108;
  sig(4,1) = 0.0005;
  sig(4,2) = 0.0018;
  sig(4,3) = 0.0035;
  sig(4,4) = 0.0034;
else;
  theta    = 0.359;
  delta    = 0.023;
  beta     = 0.99;
  gamma    = 2.33;
  pi       = 0;
  rho      = 1e-10*eye(4);
  rho(1,1) = 0.9739;
  rho0     = [1.0*(1-rho(1,1)),.195*(1-rho(2,2)),.5*(1-rho(3,3)), ...
              .23*(1-rho(4,4))]';
  sig(1,1) = 0.0096;
end;

disp('TITLE:     Tables, divisible labor cases ')
disp(' ')
help tablea
disp(' ')
disp('Parameters              Description ')
disp('--------------------------------------------------------------')
disp(sprintf('   theta    = %7.4g    capital share in production',theta))
disp(sprintf('   delta    = %7.4g    rate of capital depreciation',delta))
disp(sprintf('   beta     = %7.4g    discount factor',beta))
disp(sprintf('   gamma    = %7.4g    u(C,l)=ln(C)+gamma*ln(l)',gamma))
disp(sprintf('   pi       = %7.4g    C = c+pi*g',pi))
disp(sprintf('   rho0(1)  = %7.4g    z''=rho0+rho*z+eps,z=[lam,g,tau,phi]', ...
                                     rho0(1)))
disp(sprintf('   rho0(2)  = %7.4g    z''=rho0+rho*z+eps',rho0(2)))
disp(sprintf('   rho0(3)  = %7.4g    z''=rho0+rho*z+eps',rho0(3)))
disp(sprintf('   rho0(4)  = %7.4g    z''=rho0+rho*z+eps',rho0(4)))
disp(sprintf('   rho(1,1) = %7.4g    z''=rho0+rho*z+eps',rho(1,1)))
disp(sprintf('   rho(2,1) = %7.4g    z''=rho0+rho*z+eps',rho(1,2)))
disp(sprintf('   rho(3,1) = %7.4g    z''=rho0+rho*z+eps',rho(1,3)))
disp(sprintf('   rho(4,1) = %7.4g    z''=rho0+rho*z+eps',rho(1,4)))
disp(sprintf('   rho(2,1) = %7.4g    z''=rho0+rho*z+eps',rho(2,1)))
disp(sprintf('   rho(2,2) = %7.4g    z''=rho0+rho*z+eps',rho(2,2)))
disp(sprintf('   rho(2,3) = %7.4g    z''=rho0+rho*z+eps',rho(2,3)))
disp(sprintf('   rho(2,4) = %7.4g    z''=rho0+rho*z+eps',rho(2,4)))
disp(sprintf('   rho(3,1) = %7.4g    z''=rho0+rho*z+eps',rho(3,1)))
disp(sprintf('   rho(3,2) = %7.4g    z''=rho0+rho*z+eps',rho(3,2)))
disp(sprintf('   rho(3,3) = %7.4g    z''=rho0+rho*z+eps',rho(3,3)))
disp(sprintf('   rho(3,4) = %7.4g    z''=rho0+rho*z+eps',rho(3,4)))
disp(sprintf('   rho(4,1) = %7.4g    z''=rho0+rho*z+eps',rho(4,1)))
disp(sprintf('   rho(4,2) = %7.4g    z''=rho0+rho*z+eps',rho(4,2)))
disp(sprintf('   rho(4,3) = %7.4g    z''=rho0+rho*z+eps',rho(4,3)))
disp(sprintf('   rho(4,4) = %7.4g    z''=rho0+rho*z+eps',rho(4,4)))
disp(sprintf('   sig(1,1) = %7.4g    z''=rho0+rho*z+eps,std(eps)=sig',sig(1,1)))
disp(sprintf('   sig(2,1) = %7.4g    z''=rho0+rho*z+eps',sig(1,2)))
disp(sprintf('   sig(3,1) = %7.4g    z''=rho0+rho*z+eps',sig(1,3)))
disp(sprintf('   sig(4,1) = %7.4g    z''=rho0+rho*z+eps',sig(1,4)))
disp(sprintf('   sig(2,1) = %7.4g    z''=rho0+rho*z+eps',sig(2,1)))
disp(sprintf('   sig(2,2) = %7.4g    z''=rho0+rho*z+eps',sig(2,2)))
disp(sprintf('   sig(2,3) = %7.4g    z''=rho0+rho*z+eps',sig(2,3)))
disp(sprintf('   sig(2,4) = %7.4g    z''=rho0+rho*z+eps',sig(2,4)))
disp(sprintf('   sig(3,1) = %7.4g    z''=rho0+rho*z+eps',sig(3,1)))
disp(sprintf('   sig(3,2) = %7.4g    z''=rho0+rho*z+eps',sig(3,2)))
disp(sprintf('   sig(3,3) = %7.4g    z''=rho0+rho*z+eps',sig(3,3)))
disp(sprintf('   sig(3,4) = %7.4g    z''=rho0+rho*z+eps',sig(3,4)))
disp(sprintf('   sig(4,1) = %7.4g    z''=rho0+rho*z+eps',sig(4,1)))
disp(sprintf('   sig(4,2) = %7.4g    z''=rho0+rho*z+eps',sig(4,2)))
disp(sprintf('   sig(4,3) = %7.4g    z''=rho0+rho*z+eps',sig(4,3)))
disp(sprintf('   sig(4,4) = %7.4g    z''=rho0+rho*z+eps',sig(4,4)))
disp('--------------------------------------------------------------')

% Law of motion for x=[k,z,1, K,N,i,n]' 
% x1[t+1]=A11*x1[t]+A12*x2[t]+B1*u[t]+eps[t+1]

A11      = [1-delta,0,0,0,0,0; zeros(4,1),rho,rho0; 0,0,0,0,0,1];
A12      = [0,0;0,0;0,0;0,0;0,0;0,0];
B1       = [1,0;0,0;0,0;0,0;0,0;0,0];


% Equilibrium conditions:  x2[t] = Theta*x1[t]+Psi*u[t]

Theta    = [1,0,0,0,0,0; 0,0,0,0,0,0];
Psi      = [0,0; 0,1];


% Steady state of  [k,z,1,K,I,N,i,n]: 

z        = (eye(4)-rho)\rho0;
lam      = z(1);
g        = z(2);
tau      = z(3);
phi      = z(4);
r        = (1-beta*(1-delta+delta*tau))/(beta*(1-tau));
ydk      = r/theta;


% Calculate steady state values:

k0    = 8;
dk    = 0.0001;
for i = 1:10;
  kp     = k0 + dk;
  n0     = (ydk/lam)^(1/(1-theta))*k0;
  np     = (ydk/lam)^(1/(1-theta))*kp;
  y0     = lam*(k0/n0)^theta*n0;
  yp     = lam*(kp/np)^theta*np;
  f0     = gamma/(1-n0)-(1-phi)*(1-theta)*y0/n0/(y0-delta*k0-(1-pi)*g);
  fp     = gamma/(1-np)-(1-phi)*(1-theta)*yp/np/(yp-delta*kp-(1-pi)*g);
  fj     = (fp-f0)/dk;
  k0     = k0-fj\f0;
end;  
k        = k0;
n        = n0;
y        = y0;
z        = [k;z;1;k;n;delta*k;n];


% Calculate quadratic objective function

e        = zeros(10,1); e(6)=1;
par      = [theta;gamma;delta;pi];
v        = vfn_a(z,par);
[vz,vzz] = numder('vfn_a',z,par,ones(size(z))*1e-3,2);
M        = e*(v-vz'*z+.5*z'*vzz*z)*e'+.5*(e*vz'+vz*e'-e*z'*vzz-vzz*z*e'+vzz);
n1       = 6; 
n2       = 8; 
ns       = 10;
Q        = M(1:n2,1:n2); 
W        = M(1:n2,n2+1:ns); 
R        = M(n2+1:ns,n2+1:ns);

% Use subroutine for computing equilibrium    

[f,p]    = ddlqr(beta,A11,A12,B1,Q,R,W,Theta,Psi);


% Display answer to problem: u = -f*x1

disp(' ')
disp(' ')
disp('Intermediate results ... ')
disp(' ')
disp('   Decision Functions: ')
disp(sprintf('      i[t] =    %8.4g k[t]   + %8.4g lam[t] + %8.4g g[t] ', ...
                                 -f(1,1),-f(1,2),-f(1,3)));
disp(sprintf('              + %8.4g tau[t] + %8.4g phi[t] + %8.4g ', ...
                                 -f(1,4),-f(1,5),-f(1,6)));
disp(sprintf('      n[t] =    %8.4g k[t]   + %8.4g lam[t] + %8.4g g[t] ', ...
                                 -f(2,1),-f(2,2),-f(2,3)));
disp(sprintf('              + %8.4g tau[t] + %8.4g phi[t] + %8.4g ', ...
                                 -f(2,4),-f(2,5),-f(2,6)));
disp(' ')
disp('or ... ')
F = [f(1,1)-1+delta,f(1,2)*lam/k,f(1,3)*g/k,f(1,4)*tau/k,f(1,5)*phi/k,f(1,6)/k;
     f(2,1)*k/n,f(2,2)*lam/n,f(2,3)*g/n,f(2,4)*tau/n,f(2,5)*phi/n,f(2,6)/n]; 
disp(sprintf('      k[t+1]/k = %8.4g k[t]/k     + %8.4g lam[t]/lam + %8.4g g[t]/g', ...
                                 -F(1,1),-F(1,2),-F(1,3)));
disp(sprintf('               + %8.4g tau[t]/tau + %8.4g phi[t]/phi + %8.4g', ...
                                 -F(1,4),-F(1,5),-F(1,6)));
disp(sprintf('      n[t]/n   = %8.4g k[t]/k     + %8.4g lam[t]/lam + %8.4g g[t]/g', ...
                                 -F(2,1),-F(2,2),-F(2,3)));
disp(sprintf('              + %8.4g tau[t]/tau  + %8.4g phi[t]/phi + %8.4g', ...
                                 -F(2,4),-F(2,5),-F(2,6)));
disp(' ')


% Check on steady state

tem = -f*z(1:n1);
disp(' ')
disp('   Steady state          Decision function ')
disp('     controls        evaluated at steady state')
disp('   ------------      -------------------------')
for i=1:(ns-n2);
  disp(sprintf('  %10.4g               %10.4g  ',z(i+n2),tem(i)))
end;
disp(' ')

% For simulations use: Ao,C,x0, x=[k,lam,g,tau,phi,1,K]'

randn('seed',0)
x0     = z(1:n1);
u0     = z(n2+1:ns);
Ao     = A11+A12*Theta-A12*Psi*f-B1*f;
C      = zeros(n1); C(2:5,2:5)=sig;
nn     = 164;
nsim   = 100;
shk1   = randn(nn,nsim);
shk2   = randn(nn,nsim);
shk3   = randn(nn,nsim);
shk4   = randn(nn,nsim);
% Y  = [y,h,k,c,i,y/h,lam,g,tau,phi,w]
for i=1:nsim;
  s         = [zeros(nn,1),shk1(:,i),shk2(:,i),shk3(:,i),shk4(:,i),zeros(nn,1)];
  [x,jnk]   = edlsim(Ao,C,eye(n1),zeros(n1),s,x0);
  u         = -x*f';
  w         = (1-theta)*x(:,2).*(x(:,1).^theta.*u(:,2).^(1-theta))./u(:,2);
  o         = x(:,2).*(x(:,1).^theta.*u(:,2).^(1-theta));
  Y         = [o,u(:,2),x(:,1),o-u(:,1)-x(:,3),u(:,1),o./u(:,2),x(:,2:5),w];
  sy        = ~(mean(Y)>1e-5);
  tem       = Y+1e-8*ones(nn,1)*sy;
  Yd        = hptrend(log(tem),1600)-log(tem); 
  tem       = corrcoef(Yd);
  mn(:,i)   = mean(Y)';
  sdev(:,i) = std(Yd)'*100;
  ycor(:,i) = tem(:,1);
  ncor(:,i) = tem(:,2);
end;

disp(' ')
disp(' ')
disp('Final results ...')
disp(' ')
label = ['y  ';'n  ';'k  ';'c  ';'i  ';'y/n';'lam';'g  ';'tau';'phi';'w  '];
if nsim>1;
  tab = [mean(mn')',mean(sdev')',mean(ycor')',std(mn')', ...
         std(sdev')',std(ycor')'];
  disp('            Means      Standard Deviations    Correlation with y ')
  disp(' ---------------------------------------------------------------- ')
  for i = 1:length(label);
    disp(sprintf(['   ',label(i,:), ...
         '    %10.3g           %10.3g           %10.3g '], ...
         tab(i,1),tab(i,2),tab(i,3)))
    disp(sprintf(  ...       
         '            (%7.2g)            (%7.2g)            (%7.2g)',  ...
         tab(i,4),tab(i,5),tab(i,6)))
  end;
  disp(' ---------------------------------------------------------------- ')
  disp(' (HP filter used to calculate last 2 columns) ')
else
  tab = [mn,sdev,ycor];
  disp('            Means      Standard Deviations    Correlation with y ')
  disp(' ---------------------------------------------------------------- ')
  for i = 1:length(label);
    disp(sprintf(['   ',label(i,:), ...
         '    %10.3g           %10.3g           %10.3g '], ...
         tab(i,1),tab(i,2),tab(i,3)))
  end;
  disp(' ---------------------------------------------------------------- ')
  disp(' (HP filter used to calculate last 2 columns) ')
end;
disp(' ')
disp(' ')
tem1 = mean(ncor')';
tem2 = std(ncor')';
disp(sprintf([' Correlation between hours and productivity =', ...
     ' %10.3g (%7.2g)'],tem1(11),tem2(11)))
disp(' ')
disp(' ')
