%LSTICKY1   computes equilibria for the models in Chari, Kehoe, and 
%           McGrattan, ``Sticky Price Models of the Business Cycle: 
%           Can the Contract Multiplier Solve the Persistence Problem?''
%           This code should produce the same results as sticky1.
%           All of the first order conditions have been linearized
%           in lsticky1; only some have been linearized in sticky1.  
%
%           The optimization problems being solved are as follows.
%           Consumers solve:
%
%               max      sum sum beta^t pi(s^t) [U(C(s^t),L(s^t),M(s^t))]
%               C(s^t),   t  s^t                                 ------
%               L(s^t),                                          Pb(s^t)
%             M(s^t),B(s^t)
%
%               subj. to
%                     Pb(s^t)*C(s^t) + M(s^t)
%                          + sum_{s_{t+1}} Q(s^{t+1}|s^t)*B(s^{t+1})
%                       <= Pb(s^t)*W(s^t)*L(s^t) + M(s^{t-1}) + B(s^t) 
%                               +Profits(s^t) - Tax(s^t)
%
%           Final goods producers solve:
%
%               max      Pb - int P_i g(Y_i/Y) di
%              {Y_i}
%
%
%           Monopolists who are setting prices (group i) solve:
%
%               max   sum  Q(s^t|.)*[p*Y(i,s^t)-Pb(s^t)W(s^t)L(i,s^t)
%                p   t,s^t                     -Pb(s^t)X(i,s^t)]
%
%                    s.t.  Y(i,s^t) = D(p,Pb(s^t),{Y_j(s^t)/Y(s^t)})Y(s^t)
%                          Y(i,s^t) = F(K(i,s^t),L(i,s^t))
%                          K(i,s^t) = (1-delta) K(i,s^{t-1}) + X(i,s^t)
%                                    - phi(X(i,s^t)/K(i,s^{t-1}))*K(i,s^{t-1})
%

%           Ellen McGrattan, 2-16-99
%           Revised, ERM, 2-23-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)
%
%
%______________________________________________________________________________
%

econ     = 1;
T        = 1;
N        = 4;
psi      = 1.5;
sigma    = 1.0001;
omega    = 0.94;
eta      = 0.39;
beta     = .97^(T/N);
alp1     = 0.00001;
alp1     = 1/3;
alp2     = 1-alp1;
delta    = 1-.92^(T/N);
b        = 3;
epsilon  = 10;
chi      = epsilon+1;
rho      = 0;
rho      = 0.57^(4*T/N);
sige     = 0.00048;


% Parameters for Convex-demand economy:
if econ==2;
  chi    = -289;
end

% Parameters for Specific-factors economy:
if econ==3;
  alp1   = 0.22;
  alp2   = 0.44;
end

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-5,1e-8);
disp(sprintf(['Calculating ',int2str(nZZ),' numerical derivatives ...'],nZZ))
for i=1:nZZ;
  ZZp     = ZZ;
  ZZm     = ZZ;
  ZZp(i)  = ZZp(i)+del(i);
  ZZm(i)  = ZZm(i)-del(i);
  dR(:,i) = (res_lsticky1(ZZp,param)-res_lsticky1(ZZm,param))/(2*del(i));
end;
disp('Computing equilibrium ...')
%______________________________________________________________________________


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

disp('Generating impulse responses ...')
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)
nimp      = N-1;
npchgs    = 3/T+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]);
muprod    = 1;
pvec      = Pss*ones(1,N);

for t=1:N*npchgs;
  PMt(t,1)    = exp(G(1,:)*Xt(t,:)');                 % = P(s^t-1)/M(s^t-1)
  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}), ...
  Yt(t,1)     = exp(G(N+2,:)*Xt(t,:)');               % = Y(s^t)
  Mut(t,1)    = exp(Xt(t,nZ+1));                      % = Mu(s^t)
  Pt(t,1)     = muprod*PMt(t,1);                      % = P(s^t-1)  
                                                      %   (with M(s^0)=1)
  muprod      = muprod*Mut(t,1);                      % = Mu(s^t)..*Mu(s^1)
  Mt(t,1)     = muprod;                               % = M(s^t);
  pvec        = [pvec(2:N),Pt(t)];
  Pbt(t,1)    = gmean(pvec);                          % = Pb(s^t);
  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)
  Ct(t,1)     = Yt(t)-It(t);                          % = C(s^t)
  Yit(t,1:N)  = exp( log(Yt(t))+epsilon* ...
                (log(Pbt(t))-log(pvec(N:-1:1))) );    % = Y_1(s^t), ...
  Lit(t,1:N)  = exp( (log(Yit(t,:))- ...
                alp1*log(Kitl(t,[N,1:N-1])))/alp2 );  % = L_1(s^t), ....
  Lt(t,1)     = mean(Lit(t,1:N));                   
  MPt(t,1)    = Mt(t)/Pbt(t);                         % = M(s^t)/Pb(s^t-1);
end

figure(1)
legend off
plot(1:N*npchgs,(Yt/Yss-1)*100,'-', ...
     1:N*npchgs,(Ct/Css-1)*100,'-.', ...
     1:N*npchgs,(Lt/Lss-1)*100,'--', ...
     1:N*npchgs,(Yt/Yss-1)*100,'o', ...
     1:N*npchgs,(Ct/Css-1)*100,'x', ...
     1:N*npchgs,(Lt/Lss-1)*100,'+');
  grid
  v = axis;
  axis([N,3*N/T+N,v(3:4)])
  set(gca,'XTick',[N:N/T:3*N/T+N],'XTickLabel',[0:3])
  legend(' Output',' Consumption',' Employment')
  xlabel(['Years, (# of staggered contracts = ',int2str(N), ...
          ' and Prices Fixed ',num2str(T),' Year(s))'])
  ylabel('percent deviations')
  title('Responses to monetary disturbance')

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

disp(sprintf('The contract multiplier is equal to %6.2f ', mult));

