%TABLE2B    Use the code in sticky2.m to compute results for the 
%           near-perfect-substitutes economy in TABLE 2 of Chari, Kehoe, 
%           and McGrattan, ``Sticky Price Models of the Business 
%           Cycle: Can the Contract Multiplier Solve the Persistence 
%           Problem?''

%           Ellen McGrattan, 10-14-98
%           Revised, ERM, 2-16-99

%______________________________________________________________________________
%
%PARAMETERS:
%
%  Staggered contracts:
%
%    T = fraction of year prices are fixed (T=1 ==> fixed for 1 year)
%    N = number of cohorts
%
%  Preferences:
%
%    U(C,L,M) = [C^(1-kappa)/(1-kappa) -ome1 * L^(1+xi)/(1+xi)
%                    -ome2 * M^(1-nu)/(1-nu)]^(1-sigma)/(1-sigma)
%    beta     = the discount factor 
%    bbeta    = the discount factor used in the dynamic money demand equation
%   
%  Technology:
%
%    F(K,L)   = K^alp1 * L^alp2
%    g(Y_i/Y):  epsilon = -D'(g'(1))*g'(1), where D={g'}^(-1)
%               chi     = -D''(g'(1))*g'(1)/D'(g'(1))
%
%  Capital accumulation
%
%    K[i,t]   = (1-delta) K[i,t-1] + X[i,t] - phi(X[i,t]/K[i,t-1]) K[i,t-1]
%    phi(X/K) = b/2 * (X/K-delta)^2
%
%  Process for money:
%
%    log mu[t+1] = rho*log mu[t] + eps[t+1],   eps ~ N(0,sige)
%
%
%______________________________________________________________________________
%

%
%CASES:
%
sm       = .00001;
Params   = [ 2,.0002, sm,1-sm,          0,        sm,    .9,  1;
             2,1/.39, sm,1-sm,          0, .97^(1/8),    .9,  1;
             2,.0002,1/3, 2/3,          0,        sm,  2.76, sm;
             2,1/.39,1/3, 2/3,          0, .97^(1/8),  2.76,.028;
            13,.0002, sm,1-sm, .57^(1/13),        sm,    .9,  1;
            13,1/.39, sm,1-sm, .57^(1/13),.97^(1/52),    .9,  1;
            13,.0002,1/3, 2/3, .57^(1/13),        sm,7.0332,.0068;
            13,1/.39,1/3, 2/3, .57^(1/13),.97^(1/52),7.0332, 7];

T        = 1/4;
sigma    = 1.0001;
kappa    = .0002;
xi       = .00011;
epsilon  = 10;
chi      = 11;
sige     = 0.0001;
[nc,tem] = size(Params);

for icase=1:nc;

param    = Params(icase,:);
N        = Params(icase,1);
nu       = Params(icase,2);
alp1     = Params(icase,3);
alp2     = Params(icase,4);
rho      = Params(icase,5);
bbeta    = Params(icase,6);
ome1     = Params(icase,7);
b        = Params(icase,8);

beta     = .97^(T/N);
delta    = 1-.92^(T/N);
ome2     = 1-bbeta;

param    = [N;kappa;xi;nu;ome1;ome2;sigma;beta;bbeta; 
            alp1;alp2;delta;b;epsilon;chi];
%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

delt1    = 1-delta;
sigm1    = 1-sigma;
theta    = 1-1/epsilon;

Fk       = (1-beta*delt1)/(beta*theta);
kcratio  = 1/(Fk/alp1-delta);
lccoef   = (theta*alp2*Fk/alp1/ome1 * kcratio)^(1/(1+xi));
Css      = ((1+delta*kcratio)/(kcratio^alp1*lccoef^alp2) )^ ...
           (1/(alp1+alp2*(1-kappa)/(1+xi)-1));
Lss      = lccoef * Css^((1-kappa)/(1+xi));
Mss      = ((1-bbeta)/ome2 *Css^(-kappa))^(1/(-nu));
Kss      = (Lss^alp2*alp1/Fk)^(1/(1-alp1));
Pss      = 1/Mss;
Yss      = Lss^alp2*Kss^alp1;
%______________________________________________________________________________


%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% The ordering of variables to the residuals is
%
%     ZZ  = [z_{t+N-1},z_{t+N-2},...z_t,p_{t-2}-m_{t-2},...
%            p_{t-N}-m_{t-N},k_{1t-1},...,k_{Nt-1},mu_{t+N-1},...,mu_{t-N_1}]
%
%     z_t = [p_{t-1}-m_{t-1}
%                k_{1t}
%                  :
%                k_{Nt}
%                y_t   ]      (N+2 x 1)
%

z        = [Pss;Kss*ones(N,1);Yss];
nz       = N+2;
tem      = z*ones(1,N);
ZZ       = log([tem(:);Pss*ones(N-1,1);Kss*ones(N,1);ones(2*N-1,1)]);
nZZ      = N*N+6*N-2;
del      = max(abs(ZZ)*1e-4,1e-8);
dR       = zeros(nz,nZZ);
for i=1:nZZ;
  ZZp     = ZZ;
  ZZm     = ZZ;
  ZZp(i)  = ZZp(i)+del(i);
  ZZm(i)  = ZZm(i)-del(i);
  dR(:,i) = (res_sticky2(ZZp,param)-res_sticky2(ZZm,param))/(2*del(i));
end;
%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF DETERMINISTIC PART OF PROBLEM (A)
% 

%
% The residuals are given by
%
%     A1*[X_{t+1},Z_{t+N-1}]' + A2*[X_t,Z_{t+N-2}]' + shock terms = 0
%
% where
%
%     Z_t = [z_t,z_{t-1}, ....z_{t-N+2}]'    ((N+2)*(N-1) x 1)
%   
%     z_t = [p_{t-1}-m_{t-1}
%                k_{1t}
%                  :
%                k_{Nt}
%                y_t   ]      (N+2 x 1)
%
%     X_t = [p_{t-2}-m_{t-2}
%                  :
%            p_{t-N}-m_{t-N}
%                k_{1,t-1}
%                  :
%                k_{N,t-1}]   (2*N-1 x 1)
%
% The solution to the deterministic part of the problem is given by
%
%     Z_t = A Z_{t-1}
%
%


nX       = 2*N-1;
nZ       = nz*(N-1);
nZm      = nZ-nz;
nZp      = nZ+nz;
nXZ      = nX+nZ;
nXm      = nX+nZm;
nXp      = nX+nZp;
A1       = [eye(nX),          zeros(nX,nZ);
            zeros(nz,nX),     dR(:,1:nZ);
            zeros(nZm,nX+nz), eye(nZm)];

A2       = [zeros(1,nXm), -1,  zeros(1,nz-1);
            -eye(N-2),         zeros(N-2,nXZ-(N-2));
            zeros(N,nXZ-N-1), -eye(N),        zeros(N,1);
            dR(:,nZp+1:nXp),   zeros(nz,nZm), dR(:,nZ+1:nZp);
            zeros(nZm,nX),    -eye(nZm),      zeros(nZm,nz)];

[V,L]    = eig(A2,-A1);
L        = diag(L);
[D,I]    = sort(abs(L));
if (~( (D(nX) <= 1) & (D(nX+1) > 1) ) )
  i      = sum(D<=1);
  disp(sprintf( ...
  'WARNING: program assumes %g eigenvalues inside unit circle, found %g',nX,i));end;
V11      = V(1:nX,I(1:nX));
V21      = V(nX+1:nXZ,I(1:nX));
W        = real([V11*diag(L(I(1:nX)))/V11;
                 V21/V11]);
A        = zeros(nZ);
ii       = nXm+1:nXZ;
A(1:nz,1:N+1) ...
         = W(ii,[1,N:2*N-1]);
A(1:nz,nz+1:nz:nZ)  ...
         = W(ii,2:N-1);
A(nz+1:nZ,1:nZm) = eye(nZm);
            
%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF COEFFICIENTS ON MONETARY GROWTH (B)
%

% The residuals can be written as:
%
%   a_0*Z_{t+N-1} + a_1*Z_{t+N-2} .... + a_{N-1}*Z_t + a_N*Z_{t-1}
% + b_0*S_{t+N-1} + b_1*S_{t+N-2} .... + b_{N-1}*S_t               = 0
%
% where 
%
%     S_t = [mu_t
%            mu_{t-1}       
%                  :
%            mu_{t-N+1}]  (N x 1)
%

a0        = zeros(nz*nZ,N+1);
tem       = dR(:,1:nZ);
a0(:,1)   = tem(:);
tem       = [dR(:,nZ+1:nZp),zeros(nz,nZm)];
a0(:,N)   = tem(:);
tem       = [dR(:,[nZp+1,nZp+N:nXp]),zeros(nz,nZ-N-1)];
tem(:,nz+1:nz:nZ) ...
          = dR(:,nZp+2:nZp+N-1);
a0(:,N+1) = tem(:);

b0        = zeros(nz*N,N);
tem       = dR(:,nXp+1:nXp+N);
b0(:,1)   = tem(:);
tem       = [zeros(nz,1),dR(:,nXp+N+1:nZZ)];
b0(:,N)   = tem(:);
S         = [eye(nz);zeros(nZm,nz)];
P         = [rho,zeros(1,N-1);
             eye(N-1),zeros(N-1,1)];
D         = [rho,1,zeros(1,nZp-2);
             zeros(nZp-2,2),eye(nZp-2)];
Q         = zeros(nz,nZ);
Q(:)      = a0(:,1);
sum1      = 0;
sum2      = 0;
aa        = zeros(nz,nZ);
bb        = zeros(nz,N);
for i=1:N;
  sum1    = sum1 + kron(Q*S,(P^(N-i))');
  aa(:)   = a0(:,i+1);
  Q       = Q*A + aa;
  bb(:)   = b0(:,i);
  sum2    = sum2 + (P^(N-i))'*bb';
end;
q         = D*sum1;
d         = D*vec(sum2);
B         = -q(:,[2:nZp])\d;
B         = [0;B(1:nZp-1)];
tem       = zeros(N,nz);
tem(:)    = B;
B1        = tem';
B         = S*B1;
%
% FINAL SOLUTION:  Z[t] = A* Z[t-1] + B* S[t]
%
%______________________________________________________________________________


%______________________________________________________________________________
% 
% IMPULSE RESPONSES FOR SYSTEM
%
%       [ Z[t]   ] = [ A  B Means]  [ Z[t-1] ]      [0    ] 
%       [ S[t+1] ]   [ 0  P Means]  [ S[t]   ]   +  [Sigma] * e[t+1]
%       [   1    ]   [ 0  0   1  ]  [   1    ]      [0    ]
%
%                      (G)                      (H)         

tem       = z*ones(1,N-1);
Xbar      = log([tem(:);ones(N,1)]);
G         = [A, B; 
             zeros(N,nZ), P];
G         = [G, (eye(nZ+N)-G)*Xbar;
             zeros(1,nZ+N),1];
H         = [zeros(nZ,1); sige; zeros(N-1,1);0];

randn('seed',0)
Kit       = 0;
Kitl      = 0;
Iit       = 0;
It        = 0;
Yt        = 0;
nimp      = N-1;
npchgs    = 50/T+1;
mshk      = randn(N*npchgs,1);
imp       = [zeros(nimp-1,1);1;zeros(N*npchgs-nimp,1)];
[tem,Xt]  = dlsim(G, H, zeros(1,nZ+N+1), 0, imp, [Xbar;1]);

for t=1:N*npchgs;
  Kit(t,1:N)  = exp(G(2:N+1,:)*Xt(t,:)')';            % = K_1(s^t) ...
  Kitl(t,1:N) = exp(Xt(t,2:N+1));                     % = K_1(s^{t-1}) ...
  KK          = Kit(t,:)./Kitl(t,[N,1:N-1])-1+delta;
  Iit(t,1:N)  = (1+b*delta-sqrt(1+2*b*delta- ...
                 2*b*KK))/b .* Kitl(t,[N,1:N-1]);     % = I_1(s^t), ...
  It(t,1)     = mean(Iit(t,1:N));                     % = I(s^t)
  Yt(t,1)     = exp(G(N+2,:)*Xt(t,:)');               % = X(s^t)
end

tem       = (Yt(:)/Yt(1)-1)/(Yt(nimp+1)/Yt(1)-1);
ii        = find(tem>.5);
ll        = length(ii);
tem2      = 1-(.5-tem(ii(ll)+1))/(tem(ii(ll))-tem(ii(ll)+1));
mult      = 2*(ii(ll)-ii(1)+tem2)/N;
idivy     = (It(nimp+1)/It(1)-1)/(Yt(nimp+1)/Yt(1)-1);

Kit       = 0;
Kitl      = 0;
Iit       = 0;
It        = 0;
Yt        = 0;
In        = zeros(140,25);
Yn        = zeros(140,25);
for n=1:25;

  mshk      = randn(N*140,1);
  [tem,Xt]  = dlsim(G, H, zeros(1,nZ+N+1), 0, mshk, [Xbar;1]);

  for t=1:N*140;
    Kit(t,1:N)  = exp(G(2:N+1,:)*Xt(t,:)')';            % = K_1(s^t) ...
    Kitl(t,1:N) = exp(Xt(t,2:N+1));                     % = K_1(s^{t-1}) ...
    KK          = Kit(t,:)./Kitl(t,[N,1:N-1])-1+delta;
    Iit(t,1:N)  = (1+b*delta-sqrt(1+2*b*delta- ...
                   2*b*KK))/b .* Kitl(t,[N,1:N-1]);     % = I_1(s^t), ...
    It(t,1)     = mean(Iit(t,1:N));                     % = I(s^t)
    Yt(t,1)     = exp(G(N+2,:)*Xt(t,:)');               % = X(s^t)
  end

  II        = zeros(N,140);
  II(:)     = It;
  II        = mean(II)';
  if any(II<0 & alp1>2*sm);
    disp(['Investment below 0 on simulation ',int2str(n)])
  end;
  In(:,n)   = log(II)-hptrend(log(II),1600);
  YY        = zeros(N,140);
  YY(:)     = Yt;
  YY        = mean(YY)';
  Yn(:,n)   = log(YY)-hptrend(log(YY),1600);
end;

table2b(icase,1) = mult;
table2b(icase,2) = mean(std(In)./std(Yn));
if (alp1<2*sm); table2b(icase,2) = 0; end;
disp(sprintf('Case %2g''s multiplier and std(i)/std(y): %6.2f %6.2f', ...
     [icase,table2b(icase,:)]));

end;

table2      = zeros(4,2);
table2(:)   = table2b(:,1);
table2      = table2';
rlabel1     = 'Economy';
rlabel2     = ['Near-Perfect-Sub.';
               '                 '];
clabel3     = ['No links    ';
               'Only Money  ';
               'Only Capital';
               'Both links  '];

tab = mattab1(table2,2,0,0,0,clabel3,rlabel1,rlabel2,0);
disp('   ')
disp('   ')
disp('   ')
disp('TABLE 2b. Intertemporal Links and the Contract Multiplier');
disp(tab)

