%FIGURES    Use the code in sticky.m to compute results in the figures
%           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, 12-16-98

%______________________________________________________________________________
%
%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)
%
%
%______________________________________________________________________________
%


%
%CASES:
%
sm       = .00001;
Params   = [1, .50,    1/3,1.0001,   10,      11,     2/3, 59.1;
            1, .39,    1/3,1.0001,   10,      11,     2/3, 59.1;
            1, .35,    1/3,1.0001,   10,      11,     2/3, 59.1;
            1, .30,    1/3,1.0001,   10,      11,     2/3, 59.1;
            1, .25,    1/3,1.0001,   10,      11,     2/3, 59.1; %5
            1, .20,    1/3,1.0001,   10,      11,     2/3, 59.1;
            1, .15,    1/3,1.0001,   10,      11,     2/3, 58.8;
            1, .10,    1/3,1.0001,   10,      11,     2/3, 58.6;
            1, .05,    1/3,1.0001,   10,      11,     2/3, 58.6;
            1, .03,    1/3,1.0001,   10,      11,     2/3, 58.6; %10
            1, .01,    1/3,1.0001,   10,      11,     2/3, 56.0;
            2, .39,   .001,1.0001,   10,      11,    .999,193.3;
            2, .39,    .01,1.0001,   10,      11,     .99,191.3; 
            2, .39,    .05,1.0001,   10,      11,     .95,177.0;
            2, .39,    .10,1.0001,   10,      11,     .90,157.5; %15
            2, .39,    .15,1.0001,   10,      11,     .85,137.5;
            2, .39,    .20,1.0001,   10,      11,     .80,116.5;
            2, .39,    .25,1.0001,   10,      11,     .75, 95.5; 
            2, .39,    .30,1.0001,   10,      11,     .70, 73.9;
            2, .39,    1/3,1.0001,   10,      11,     2/3, 59.1; %20
            2, .39,    .35,1.0001,   10,      11,     .65, 51.2;
            2, .39,    .40,1.0001,   10,      11,     .60, 27.5;
            3, .39,    1/3,   .01,   10,      11,     2/3,111.0; 
            3, .39,    1/3,    .1,   10,      11,     2/3,106.5;
            3, .39,    1/3,    .2,   10,      11,     2/3,101.0; %25
            3, .39,    1/3,    .3,   10,      11,     2/3, 95.5;
            3, .39,    1/3,    .5,   10,      11,     2/3, 85.0;
            3, .39,    1/3,   .75,   10,      11,     2/3, 72.0; 
            3, .39,    1/3,1.0001,   10,      11,     2/3, 59.1;
            3, .39,    1/3,  1.25,   10,      11,     2/3, 46.0; %30
            3, .39,    1/3,   1.5,   10,      11,     2/3, 32.5;
            3, .39,    1/3,  1.75,   10,      11,     2/3, 20.0;
            3, .39,    1/3,  2.00,   10,      11,     2/3, 10.0; 
            3, .39,    1/3,  2.25,   10,      11,     2/3,  5.7;
            3, .39,    1/3,   2.7,   10,      11,     2/3,  3.5; %35
            3, .39,    1/3,  2.75,   10,      11,     2/3,  3.5;
            3, .39,    1/3,     3,   10,      11,     2/3,  3.0;
            3, .39,    1/3,     4,   10,      11,     2/3,  3.0;
            3, .39,    1/3,     5,   10,      11,     2/3,  3.0;
            3, .39,    1/3,    10,   10,      11,     2/3,  3.0; %40
            4, .39,    1/3,1.0001,    5,       6,     2/3, 76.0;
            4, .39,    1/3,1.0001,    6,       7,     2/3, 70.5;
            4, .39,    1/3,1.0001,   10,      11,     2/3, 59.1; 
            4, .39,    1/3,1.0001,   20,      21,     2/3, 49.4;
            4, .39,    1/3,1.0001,   30,      31,     2/3, 45.4; %45
            4, .39,    1/3,1.0001,   40,      41,     2/3, 43.5;
            4, .39,    1/3,1.0001,   50,      51,     2/3, 41.8;
            4, .39,    1/3,1.0001,  100,     101,     2/3, 37.5; 
            4, .39,    1/3,1.0001,  150,     151,     2/3  35.5;
            4, .39,    1/3,1.0001,  200,     201,     2/3, 34.0; %50
            4, .39,    1/3,1.0001,  250,     251,     2/3,   33;
            4, .39,    1/3,1.0001,  300,     301,     2/3, 32.3;
            5, .39,    1/3,1.0001,   10,      11,     2/3, 59.1; 
            5, .39,    1/3,1.0001,   10,    -289,     2/3, 52.0;
            5, .39,    1/3,1.0001,   10,  -289*2,     2/3, 49.1; %55
            5, .39,    1/3,1.0001,   10,  -289*3,     2/3, 47.5;
            5, .39,    1/3,1.0001,   10,  -289*4,     2/3, 46.0;
            5, .39,    1/3,1.0001,   10,  -289*5,     2/3, 45.0; 
            5, .39,    1/3,1.0001,   10, -289*10,     2/3, 41.5;
            5, .39,    1/3,1.0001,   10, -289*50,     2/3, 34.8; %60
            5, .39,    1/3,1.0001,   10,-289*100,     2/3, 33.0;
            5, .39,    1/3,1.0001,   10,-289*150,     2/3, 32.1;
            5, .39,    1/3,1.0001,   10,-289*200,     2/3, 31.9; 
            6, .39,    1/3,1.0001,   10,      11,     2/3, 59.1; 
            6, .39,1/3/1.1,1.0001,   10,      11, 2/3/1.1, 72.0; %65
            6, .39,1/3/1.2,1.0001,   10,      11, 2/3/1.2, 82.5; 
            6, .39,1/3/1.3,1.0001,   10,      11, 2/3/1.3,   91;
            6, .39,1/3/1.5,1.0001,   10,      11, 2/3/1.5,105.5; 
            6, .39,1/3/1.7,1.0001,   10,      11, 2/3/1.7,116.5;
            6, .39,1/3/2.0,1.0001,   10,      11, 2/3/2.0,  130; %70
            7, .39,    1/3,1.0001,2.667,    -289,     2/3,   80; 
            7, .39,    1/3,1.0001,    4,    -289,     2/3, 67.5;
            7, .39,    1/3,1.0001,    5,    -289,     2/3, 62.5; 
            7, .39,    1/3,1.0001,    6,    -289,     2/3, 59.1;
            7, .39,    1/3,1.0001,   10,    -289,     2/3, 52.1; %75
            7, .39,    1/3,1.0001,   20,    -289,     2/3, 46.1; 
            7, .39,    1/3,1.0001,   30,    -289,     2/3, 43.5;
            7, .39,    1/3,1.0001,   40,    -289,     2/3,   42; 
            7, .39,    1/3,1.0001,   50,    -289,     2/3,   41;
            7, .39,    1/3,1.0001,  100,    -289,     2/3, 37.5; %80
            7, .39,    1/3,1.0001,  150,    -289,     2/3  35.5; 
            7, .39,    1/3,1.0001,  200,    -289,     2/3,   34;
            7, .39,    1/3,1.0001,  250,    -289,     2/3,   33; 
            7, .39,    1/3,1.0001,  300,    -289,     2/3, 32.5;
            8, .39,    .22,1.0001,2.667,   3.667,     .44,135.5; %85
            8, .39,    .22,1.0001,    4,       5,     .44,  123; 
            8, .39,    .22,1.0001,    5,       6,     .44,  118;
            8, .39,    .22,1.0001,    6,       7,     .44,  114; 
            8, .39,    .22,1.0001,   10,      11,     .44,106.5;
            8, .39,    .22,1.0001,   20,      21,     .44, 99.5; %90 
            8, .39,    .22,1.0001,   30,      31,     .44, 96.5; 
            8, .39,    .22,1.0001,   40,      41,     .44, 94.5;
            8, .39,    .22,1.0001,   50,      51,     .44, 93.0; 
            8, .39,    .22,1.0001,  100,     101,     .44, 88.8;
            8, .39,    .22,1.0001,  150,     151,     .44  86.5; %95
            8, .39,    .22,1.0001,  200,     201,     .44, 86.0;
            8, .39,    .22,1.0001,  250,     251,     .44, 86.0;
            8, .39,    .22,1.0001,  300,     301,     .44, 87.5]; 

N        = 13;
T        = 1/4;
psi      = 1.5;
omega    = 0.94;
rho      = 0.57^(4*T/N);
sige     = 0.0001;
[nc,tem] = size(Params);

for icase=1:nc;

econ     = Params(icase,1);
eta      = Params(icase,2);
alp1     = Params(icase,3);
sigma    = Params(icase,4);
epsilon  = Params(icase,5);
chi      = Params(icase,6);
alp2     = Params(icase,7);
b        = Params(icase,8);
beta     = .97^(T/N);
delta    = 1-.92^(T/N);

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);
    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;

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

end;  

i1          = find(Params(:,1)==1);
i2          = find(Params(:,1)==2);
i3          = find(Params(:,1)==3);
i4          = find(Params(:,1)==4);
i5          = find(Params(:,1)==5);
i6          = find(Params(:,1)==6);
i7          = find(Params(:,1)==7);
i8          = find(Params(:,1)==8);

x1          = Params(i1,2);
x2          = Params(i2,3);
x3          = Params(i3,4);
x4          = Params(i4,5);
x5          = Params(i5,6);
x6          = 1-Params(i6,3)-Params(i6,7);
x7          = Params(i7,5);
x8          = Params(i8,5);

y1          = tab(i1,1);
y2          = tab(i2,1);
y3          = tab(i3,1);
y4          = tab(i4,1);
y5          = tab(i5,1);
y6          = tab(i6,1);
y7          = tab(i7,1);
y8          = tab(i8,1);

figure(1)
subplot(2,2,1), plot(x1,y1),axis([0,.5,0,2]),  ...
   title('Sensitivity analysis for the benchmark economy.'),  ...
   xlabel('A. Interest elasticity'), ylabel('Contract Multiplier')
subplot(2,2,2), plot(x2,y2),axis([0,.4,0,2]),  ...
   xlabel('B. Capital share'),       ylabel('Contract Multiplier')
subplot(2,2,3), plot(x3,y3),axis([0,10,0,2]),  ...
   xlabel('C. Risk aversion'),       ylabel('contract multiplier')
subplot(2,2,4), plot(x4,y4),axis([0,300,0,2]),  ...
   xlabel('D. Demand elasticity'),   ylabel('contract multiplier')

figure(2)
subplot(2,2,1), plot(x5,y5),axis([-50000,0,0,5]),  ...
   title('Sensitivity analysis for two economies.'), ...
   xlabel('A. Demand Curvature, Convex Demand'), ...
   ylabel('Contract Multiplier')
subplot(2,2,2), plot(x6,y6),axis([0,.5,0,2]), ...
   xlabel('B. Fixed Factor''s Share of Income, Specific Factors'), ...
   ylabel('Contract Multiplier')
subplot(2,2,3),plot(x7,y7),axis([0,300,1,4]),  ...
   xlabel('C. Demand Elasticity, Convex Demand'), ...
   ylabel('Contract Multiplier')
subplot(2,2,4),plot(x8,y8),axis([0,300,1,4]),  ...
   xlabel('D. Demand Elasticity, Specific Factors'), ...
   ylabel('Contract Multiplier')
