%TABLE3     Use the code in sticky1.m to compute results in TABLE 3
%           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) = [CES(C,M) * (1-L)^psi]^(1-sigma)/(1-sigma)
%    CES(C,M) = [omega*C^((eta-1)/eta)+(1-omega)*M^((eta-1)/eta)]^(eta/(eta-1))
%    beta     = the discount factor 
%   
%  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)
%
%
%  Benchmark parameterizations:
%
%    T       = 13
%    N       = 1/4
%    omega   = 0.94
%    eta     = 0.39
%    sigma   = 1.0001
%    psi     = 1.5
%    beta    = .97^(T/N)
%    alp1    = 1/3            (or 2/9 with specific factors)
%    alp2    = 2/3            (or 4/9 with specific factors)
%    epsilon = 10
%    chi     = 11             (or -289 with convex demand)
%    delta   = 1-.92^(T/N)
%    rho     = 0.57^(4*T/N)
%
%  Parameter search:
%
%  (A) eta     in [  .2,  .5]
%      alp1    in [ .22, .44] 
%      psi     in [ .11,   9]
%      sigma   in [   1,  10]
%      epsilon in [   5, 300]
%      rho     in [   0, .99]
%
%  (B) Add:
%      chi     in [-289,  11]
%
%  (C) Add:
%      alp1    in [ .22, .44]
%      alp2    in [ .33, .78]
%______________________________________________________________________________
%

%
%CASES:
%
sm       = .00001;
Params   = [ 1,1, .2,300,301,.22,.78,1.0001,  9,.10^(1/13),  88.5; 
             2,1, .2,300,301,.22,.78,1.0001,  9,.20^(1/13),  88.5; 
             3,1, .2,300,301,.22,.78,1.0001,  9,.30^(1/13),  88.5; 
             4,1, .2,300,301,.22,.78,1.0001,  9,.40^(1/13),  88.5; 
             5,1, .2,300,301,.22,.78,1.0001,  9,.50^(1/13),  88.5; 
             6,1, .2,300,301,.22,.78,1.0001,  9,.57^(1/13),  88.5; 
             7,1, .2,300,301,.22,.78,1.0001,  9,.60^(1/13),  88.4;
             8,1, .2,300,301,.22,.78,1.0001,  9,.70^(1/13),  88.4; 
             9,1, .2,300,301,.22,.78,1.0001,  9,.80^(1/13),  89.0;
            10,1, .2,300,301,.22,.78,1.0001,  9,.90^(1/13),  91.9;
            11,1,.25,300,301,.22,.78,1.0001,  9,.60^(1/13),  88.5; 
            12,1, .2,280,281,.22,.78,1.0001,  9,.60^(1/13),  88.7; 
            13,1, .2,300,301,.33,.67,1.0001,  9,.60^(1/13),  35.2;
            14,1, .2,300,301,.22,.78,1.5001,  9,.60^(1/13),  83.5;
            15,1, .2,300,301,.22,.78,1.0001,  8,.60^(1/13),  88.0;
            16,2  .2, 5,-289,.22,.78,1.0001,  9,.10^(1/13), 109.5;
            17,2  .2, 5,-289,.22,.78,1.0001,  9,.20^(1/13), 109.5;
            18,2  .2, 5,-289,.22,.78,1.0001,  9,.30^(1/13), 109.5;
            19,2  .2, 5,-289,.22,.78,1.0001,  9,.40^(1/13), 109.5;
            20,2  .2, 5,-289,.22,.78,1.0001,  9,.50^(1/13), 109.5; 
            21,2  .2, 5,-289,.22,.78,1.0001,  9,.60^(1/13), 110.0; 
            22,2  .2, 5,-289,.22,.78,1.0001,  9,.70^(1/13), 110.0; 
            23,2  .2, 5,-289,.22,.78,1.0001,  9,.80^(1/13), 110.4; 
            24,2  .2, 5,-289,.22,.78,1.0001,  9,.90^(1/13), 112.3; 
            25,2 .25, 5,-289,.22,.78,1.0001,  9,.70^(1/13), 110.2; 
            26,2  .2, 6,-289,.22,.78,1.0001,  9,.70^(1/13), 107.5; 
            27,2  .2, 5,-250,.22,.78,1.0001,  9,.70^(1/13), 110.7; 
            28,2  .2, 5,-289,.33,.67,1.0001,  9,.70^(1/13),  67.0; 
            29,2  .2, 5,-289,.22,.78,1.5001,  9,.70^(1/13), 107.5; 
            30,2  .2, 5,-289,.22,.78,1.0001,  8,.70^(1/13), 109.0; 
            31,3, .2,300,301,.22,.33,1.0001,  9,.10^(1/13),  88.5; 
            32,3, .2,300,301,.22,.33,1.0001,  9,.20^(1/13),  88.2; 
            33,3, .2,300,301,.22,.33,1.0001,  9,.30^(1/13),  87.9; 
            34,3, .2,300,301,.22,.33,1.0001,  9,.40^(1/13),  87.6; 
            35,3, .2,300,301,.22,.33,1.0001,  9,.50^(1/13),  87.3; 
            36,3, .2,300,301,.22,.33,1.0001,  9,.60^(1/13),  87.0; 
            37,3, .2,300,301,.22,.33,1.0001,  9,.70^(1/13),  86.4; 
            38,3, .2,300,301,.22,.33,1.0001,  9,.80^(1/13),  86.1; 
            39,3, .2,300,301,.22,.33,1.0001,  9,.90^(1/13),  87.5; 
            40,3,.25,300,301,.22,.33,1.0001,  9,.70^(1/13),  87.0; 
            41,3, .2,280,281,.22,.33,1.0001,  9,.70^(1/13),  86.5; 
            42,3, .2,  5,  6,.22,.33,1.0001,  9,.70^(1/13), 119.5; 
            43,3, .2,300,301,.22,.44,1.0001,  9,.70^(1/13),    84; 
            44,3, .2,300,301,.33,.33,1.0001,  9,.70^(1/13),  36.7; 
            45,3, .2,300,301,.22,.33,1.5001,  9,.70^(1/13),    82; 
            46,3, .2,300,301,.22,.33,1.0001,  8,.70^(1/13),  86.2; 
            47,3, .2,300,301,.22,.33,1.0001,.11,.70^(1/13),   113];


N        = 13;
T        = 1/4;
beta     = .97^(T/N);
delta    = 1-.92^(T/N);
omega    = 0.94;
sige     = 0.0001;
[nc,tem] = size(Params);

for icase=1:nc;

eta      = Params(icase,3);
epsilon  = Params(icase,4);
chi      = Params(icase,5);
alp1     = Params(icase,6);
alp2     = Params(icase,7);
sigma    = Params(icase,8);
psi      = Params(icase,9);
rho      = Params(icase,10);
b        = Params(icase,11);

param    = [N;psi;sigma;omega;eta;beta;alp1;alp2;delta;b;epsilon;chi];
%______________________________________________________________________________


%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

omeg1    = 1-omega;
delt1    = 1-delta;
sigm1    = 1-sigma;
theta    = 1-1/epsilon;
nu       = (eta-1)/eta;
Fk       = (1-beta*delt1)/(beta*theta);
ckratio  = Fk/alp1-delta;
mcratio  = (omega/omeg1*(1-beta))^(-eta);
Lss      = 1/(1+psi*ckratio*(1+omeg1/omega*mcratio^(1-1/eta))/ ...
           (theta*Fk*alp2/alp1));
Kss      = (alp1*Lss^alp2/Fk)^(1/(1-alp1));
Css      = ckratio*Kss;
Mss      = mcratio*Css;
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_sticky1(ZZp,param)-res_sticky1(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    = 3/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,:)');               % = Y(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);
    error(['Investment below 0 on simultion ',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;

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

end;  

i1       = find(Params(:,2)==1);
i2       = find(Params(:,2)==2);
i3       = find(Params(:,2)==3);
j1       = i1( find(max(table3(i1,1)) == table3(i1,1)) );
j2       = i2( find(max(table3(i2,1)) == table3(i2,1)) );
j3       = i3( find(max(table3(i3,1)) == table3(i3,1)) );

tab3     = [Params(j1,3:11),table3(j1,1);
            Params(j2,3:11),table3(j2,1);
            Params(j3,3:11),table3(j3,1)]';


clabel1  = 'Economy';
clabel3  = ['   Benchmark    ';
            ' Convex Demand  ';
            'Specific Factors'];
rlabel1  = ['          ';
            'Parameter '];
rlabel2  = ['eta        ';
            'epsilon    ';
            'chi        ';
            'alp1       ';
            'alp2       ';
            'sigma      ';
            'psi        ';
            'rho        ';
            'b          ';
            'MULTIPLIER '];

tab      = mattab1(tab3,2,0,clabel1,0,clabel3,rlabel1,rlabel2,0);
disp('   ')
disp('   ')
disp('   ')
disp('   ')
disp('TABLE 3. Maximum Multipliers Found in Three Economies')
disp(tab)

