%TRIPLE       computes equilibria for a ``triple-sticky'' business cycle 
%             model with: staggered sticky prices and wages, complete 
%             markets, habit persistence, and multiple shocks.  
%
%             The optimization problems being solved are as follows.
%             UNION j solves:
%
%               max       sum  beta^t pi(s^t) *
%              C(j,s^t)  t,s^t     U[C(j,s^t),C(j,s^{t-1}),L(j,s^t),M(j,s^t)]
%              M(j,s^t)                                             --------
%              B(j,s^t)                                              P(s^t)
%              W(j,s^{t-1})
%
%             subject  to
%                 P(s^t)*C(j,s^t) + M(j,s^t) 
%                     + sum_{s_{t+1}} Q(s^{t+1}|s^t)*B(j,s^{t+1})
%                   <=  P(s^t)*W(j,s^{t-1})*L(j,s^t)) + M(j,s^{t-1}) 
%                          + B(j,s^t) + Profits(s^t) - Tax(s^t)
%
%                 L(j,s^t)   = int l(i,j,s^t) di
%                 l(i,j,s^t) = (Wbar(s^t)/W(j,s^{t-1}))^(1/(1-v)) L^d(i,s^t)
%
%
%             FINAL GOODS PRODUCERS solve:
%
%               max     P*(int Y(i)^theta di)^{1/theta} -int P(i) Y(i) di 
%               Y(i)  
%
%             INTERMEDIATE GOODS PRODUCERS solve:
%
%               max    sum  Q(s^t|.)*[P(i,s^t)*Y(i,s^t)
%               P(i)   t,s^t           - P(s^t) int W(j,s^t)l(i,j,s^t) dj
%                                      - P(s^t) X(i,s^t)]
%
%                       s.t. Y(i,s^t)   = F(K(i,s^{t-1}),A(s^t) L^d(i,s^t))
%                            L^d(i,s^t) = [ int l(i,j,s^t)^v dj ] ^(1/v)
%
%                            Y(i)       = Y (P/P(i))^(1/(1-theta))
%                            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})
%
%                            P(i,s^tau) = P(i,s^{t-1}), tau=t,...,t+N-1  
%                            P(i,s^tau) = P(i,s^{t+N-1}), tau=t+N,...,t+2N-1  
%                                        :
%           INTEREST RATES:
%
%              r(s^t) = irule' [     r(s^{t-1})
%                                    r(s^{t-2})
%                                    r(s^{t-3})
%                                E_t log P(s^{t+1})-log P(s^t)
%                                    log P(s^t)    -log P(s^{t-1})
%                                    log P(s^{t-1})-log P(s^{t-2})
%                                    log P(s^{t-2})-log P(s^{t-3})
%                                    log y(s^t)
%                                    log y(s^{t-1})
%                                    log y(s^{t-2})  ]  + eps_r(s^t)
%
%

%             Ellen McGrattan, 3-21-00
%             Revised, ERM, 3-20-04

%______________________________________________________________________________
%
%PARAMETERS:
%
%  Staggered contracts:
%
%    T  = fraction of year prices are fixed (T=1 ==> fixed for 1 year)
%    N  = number of price-setting cohorts 
%    Nw = number of wage-setting cohorts
%
%  Preferences:
%
%    Utype=1:
%    U(C,Clag,M,L) = [CES(C,M) * (1-L)^psi]^(1-sigma)/(1-sigma)
%
%    Utype=2:
%    U(C,Clag,M,L) = CES(C,M)^(1-sigma)/(1-sigma)+psi*(1-L)^(1-xi)/(1-xi)
%
%    Utype=3:
%    U(C,Clag,M,L) = [C^(1-kappa)/(1-kappa) -psi * L^(1+xi)/(1+xi)
%                                    +omega * M^nu/nu]^(1-sigma)/(1-sigma)
%    Utype=4:
%    U(C,Clag,M,L) = [C^(1-kappa)/(1-kappa) +psi * (1-L)^(1-xi)/(1-xi)
%                                    +omega * M^nu/nu]^(1-sigma)/(1-sigma)
%    Utype=5:
%    U(C,Clag,M,L) = shock*log(C-kappa*Clag)-psi * L^(1+xi)/(1+xi)
%                                    + omega * M^nu/nu
%    where
%      CES(C,M) = [omega*C^nu+(1-omega)*M^nu]^(1/nu)
%      nu       = (eta-1)/eta
%
%  Discounting:
%
%    beta     = the discount factor
%
%  Technology:
%
%    F(K,AL)  = K^alpha * (AL)^(1-alpha)
%    Y        = (int Y(i)^theta di)^(1/theta)
%
%  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
%
%  Taylor rule for interest rates:
%
%    r[t]     = irule' * [rhs variables ] + constant + eps_r[t]
%    mu       = mean rate of inflation and money growth
%
%  Timing of shocks:
%
%    Price and wage-setters observe only lagged shocks when choosing prices
%  
%______________________________________________________________________________
%

T        = 1;
N        = 4;
Nw       = 2;
Utype    = 3;
eta      = .39;
kappa    = .42;
omega    = 1-.97^(T/N);
psi      = 1;
sigma    = 0.0000000001;
xi       = .8;
Utype    = 5;
upars    = [eta,kappa,omega,psi,sigma,xi];
beta     = .97^(T/N);
alpha    = 1/3;
delta    = 1-.92^(T/N);
b        = 20;
theta    = .9;
v        = .87;
mu       = 1.0;
gs       = 0.6;
rhor     = 0;
rhog     = .95;
rhoa     = .95;
rhop     = .0;
sigr     = 1;
sigg     = 1;
siga     = 1;
sigp     = 1;
irule    = [0.79;0;0; .21*2.15;0;0;0; .21*.93/4;0;0];
irule    = irule+1e-3;  % needed to aid computation of generalized eigs
alph1    = 1-alpha;
thet1    = 1-theta;
v1       = 1-v;
param    = [T,N,Nw,Utype,upars,beta,alpha,delta,b,theta,v,mu,gs,irule']';


%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

ss       = tripless(param);
param    = [param;ss];

%______________________________________________________________________________



%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% Let mN = max(N,Nw).
% The ordering of variables to the residuals is
%
%     ZZ  = [z_{t+mN-1},z_{t+mN-2},...z_t,X_t,
%            eps_{r,t},g_{t+mN-1},...,g_t,
%            a_{t+mN-1},...,a_t,
%            pshk_{t+mN-1},...,pshk_t]'
%
%     z_t = [p_{t-1}           NOTES: (1) variables are logged and demeaned
%            w_{t-1}                  (2) prices are normalized by mu^t
%            k_{1t}   
%              :
%            k_{Nt}
%            c_{1t}
%              :
%            c_{Nw,t}
%            md_{1t}
%              :
%            md_{Nw,t}
%            y_t
%            r_t ]      (N+2*Nw+4  x 1)
%

ps       = ss(1);
ws       = ss(2);
ks       = ss(3:N+2);
cs       = ss(N+3:N+Nw+2);
mds      = ss(N+Nw+3:N+2*Nw+2);
ys       = ss(N+2*Nw+3);
rs       = mu/beta-1;
z        = [ps;ws;ks;cs;mds;ys;exp(rs)];
X        = [ps*ones(N+2,1);ws*ones(Nw-1,1);ks;cs;ys;ys;exp(rs);exp(rs);exp(rs)];
nz       = N+2*Nw+4;
mN       = max(N,Nw);
tem      = z*ones(1,mN);
ZZ       = log([tem(:);X;exp(0);gs*ones(mN,1);ones(mN,1);ones(mN,1)]);
nZZ      = length(ZZ);
del      = max(abs(ZZ)*1e-4,1e-8);
for i=1:nZZ;
  ZZp     = ZZ;
  ZZm     = ZZ;
  ZZp(i)  = ZZp(i)+del(i);
  ZZm(i)  = ZZm(i)-del(i);
  dR(:,i) = (res_triple(ZZp,param)-res_triple(ZZm,param))/(2*del(i));
end;
%______________________________________________________________________________


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

%
% The residuals are given by
%
%     A1*[X_{t+1},z_{t+mN-1},..., z_{t+1}]'
%   + A2*[X_t,    z_{t+mN-2},..., z_t]'     + shock terms = 0
%
% where
%
%     z_t = [p_{t-1}
%            w_{t-1}
%            k_{1t}
%              :
%            k_{Nt}
%            c_{1t}
%              :
%            c_{Nw,t}
%            md_{1t}
%              :
%            md_{Nw,t}
%            y_t
%            r_t ]      (N+2*Nw+4  x 1)
%
%     X_t = [p_{t-2}
%                :
%            p_{t-N}
%            p_{t-N-1}
%            p_{t-N-2}
%            p_{t-N-3}
%            w_{t-2}
%                :
%            w_{t-Nw}
%            k_{1,t-1}
%                :
%            k_{N,t-1}
%            c_{1,t-1}
%             :
%            c_{Nw,t-1}
%            y_{t-1}
%            y_{t-2}
%            r_{t-1}
%            r_{t-2}
%            r_{t-3}]        (2*N+2*Nw+6 x 1)
%
% The solution to the deterministic part of the problem is given by
%
%     Z_t = A Z_{t-1}
%
% where
%
%     Z_t = [z_t,z_{t-1}, ....z_{t-N-1}]'    (nz*(mN+2) x 1)
%
%

nX       = 2*N+2*Nw+6;
nZ       = (mN+2)*nz;
nA1      = nX+(mN-1)*nz;
A1       = [eye(nX,nX),zeros(nX,(mN-1)*nz);
            zeros(nz,nX),dR(:,1:(mN-1)*nz);
            zeros((mN-2)*nz,nX+nz),eye((mN-2)*nz)];

I1       = [diag(ones(N+1,1),-1),zeros(N+2,nX-N-2);
            zeros(Nw-1,N+2),diag(ones(Nw-2,1),-1),zeros(Nw-1,N+Nw+5);
            zeros(N+Nw,nX);
            zeros(5,nX-5),[0,0,0,0,0;
                           1,0,0,0,0;
                           0,0,0,0,0;
                           0,0,1,0,0;
                           0,0,0,1,0]];

I2       = [1,zeros(1,nz-1);
            zeros(N+1,nz);
            0,1,zeros(1,nz-2);
            zeros(Nw-2,nz);
            zeros(N+Nw,2),eye(N+Nw),zeros(N+Nw,Nw+2);
            zeros(5,nz-2),[1,0;
                           0,0;
                           0,1;
                           0,0;
                           0,0]];

A2       = [-I1,zeros(nX,(mN-2)*nz),-I2;
            dR(:,mN*nz+1:mN*nz+nX),zeros(nz,(mN-2)*nz), ...
            dR(:,(mN-1)*nz+1:mN*nz);
            zeros((mN-2)*nz,nX),-eye((mN-2)*nz),zeros((mN-2)*nz,nz)];

[V,L]    = eig(A2,-A1);
L        = diag(L);
[D,I]    = sort(abs(L));
if (~( (D(nX) <= 1) & (D(nX+1) >= 1) ) )
  disp('Can''t order eigenvalues')
end;
V11      = V(1:nX,I(1:nX));
V21      = V(nX+1:nA1,I(1:nX));
W        = real([V11*diag(L(I(1:nX)))/V11;
                 V21/V11]);

A                         = zeros((mN+2)*nz);
A(1:nz,1:nz:(N+2)*nz)     = W(nA1-nz+1:nA1,1:N+2);
A(1:nz,2:nz:(Nw-1)*nz)    = W(nA1-nz+1:nA1,N+3:N+Nw+1);
A(1:nz,3:N+2)             = W(nA1-nz+1:nA1,N+Nw+2:2*N+Nw+1);
A(1:nz,N+3:N+Nw+2)        = W(nA1-nz+1:nA1,2*N+Nw+2:2*N+2*Nw+1);
A(1:nz,nz-1)              = W(nA1-nz+1:nA1,2*N+2*Nw+2);
A(1:nz,2*nz-1)            = W(nA1-nz+1:nA1,2*N+2*Nw+3);
A(1:nz,nz)                = W(nA1-nz+1:nA1,2*N+2*Nw+4);
A(1:nz,2*nz)              = W(nA1-nz+1:nA1,2*N+2*Nw+5);
A(1:nz,3*nz)              = W(nA1-nz+1:nA1,2*N+2*Nw+6);
A(nz+1:nZ,1:nZ-nz)        = eye(nZ-nz);

%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF COEFFICIENTS B
%

% The residuals can be written as:
%
%   a_0*Z_{t+mN-1} + a_1*Z_{t+mN-2} .... + a_{mN-1}*Z_t + a_mN*Z_{t-1}
% + b_0*S_{t+mN-1} + b_1*S_{t+mN-2} .... + b_{mN-1}*S_t               = 0
%
% where
%
%     S_t = [eps_{r,t}
%            eps_{r,t-1}
%            g_t
%            g_{t-1}
%            a_t
%            a_{t-1}
%            pshk_t
%            pshk_{t-1}]     (8 x 1)
%

nS             = 8;
a0             = zeros(nz*nZ,mN+1);
b0             = zeros(nz*nS,mN);
a0(:,1)        = vec([dR(:,1:mN*nz),zeros(nz,2*nz)]);
tem            = zeros(nz,nZ);
tem(:,1:nz:(N+2)*nz)  = dR(:,mN*nz+1:mN*nz+N+2);
tem(:,2:nz:(Nw-1)*nz) = dR(:,mN*nz+N+3:mN*nz+N+Nw+1);
tem(:,3:N+2)          = dR(:,mN*nz+N+Nw+2:mN*nz+2*N+Nw+1);
tem(:,N+3:N+Nw+2)     = dR(:,mN*nz+2*N+Nw+2:mN*nz+2*N+2*Nw+1);
tem(:,nz-1:nz:2*nz-1) = dR(:,mN*nz+2*N+2*Nw+2:mN*nz+2*N+2*Nw+3);
tem(:,nz:nz:3*nz)     = dR(:,mN*nz+2*N+2*Nw+4:mN*nz+2*N+2*Nw+6);
a0(:,mN+1)     = vec(tem);
for i=1:mN;
  b0(:,i)      = vec([zeros(nz,2),dR(:,mN*nz+nX+i+1),zeros(nz,1), ...
                                  dR(:,mN*nz+nX+mN+i+1),zeros(nz,1), ...
                                  dR(:,mN*nz+nX+2*mN+i+1),zeros(nz,1)]);
end;
b0(:,mN)       = vec([dR(:,mN*nz+nX+1),zeros(nz,7)]);
sb             = [eye(nz);zeros(nZ-nz,nz)];
sa             = zeros(nz,nZ);
sa(:)          = a0(:,1);
P              = [rhor,0,0,0,0,0,0,0;
                  1,0,0,0,0,0,0,0;
                  0,0,rhog,0,0,0,0,0;
                  0,0,1,0,0,0,0,0;
                  0,0,0,0,rhoa,0,0,0;
                  0,0,0,0,1,0,0,0;
                  0,0,0,0,0,0,rhop,0;
                  0,0,0,0,0,0,1,0];

sum1           = 0;
sum2           = 0;
aa             = zeros(nz,nZ);
bb             = zeros(nz,nS);
for i=1:mN;
  sum1         = sum1 + kron(sa*sb,(P^(mN-i))');
  aa(:)        = a0(:,i+1);
  sa           = sa*A + aa;
  bb(:)        = b0(:,i);
  sum2         = sum2 + (P^(mN-i))'*bb';
end;

Psi1           = [rhor,1,0,0,0,0,0,0;
                  0,0,rhog,1,0,0,0,0;
                  0,0,0,0,rhoa,1,0,0;
                  0,0,0,0,0,0,rhop,1];
dd             = [kron(eye(2),Psi1),zeros(8,nz*nS-16);
                  zeros(nz*nS-16,16),eye(nz*nS-16)];


nzero          = [2:2:16,17:nz*nS];
qq             = dd*sum1;
rr             = dd*vec(sum2);
tem            = zeros(nz*nS,1);
tem(nzero)     = -qq(:,nzero)\rr;
B1             = zeros(nS,nz);
B1(:)          = tem;
B1             = B1';
B              = sb*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,mN+2);
Xbar      = log([tem(:);1;1;gs;gs;1;1;1;1]);
G         = [A, B;
             zeros(nS,nZ), P];
G         = [G, (eye(nZ+nS)-G)*Xbar;
             zeros(1,nZ+nS),1];
H         = [zeros(nZ,4);
             sigr,0,0,0;
             0,0,0,0;
             0,sigg,0,0;
             0,0,0,0;
             0,0,siga,0;
             0,0,0,0;
             0,0,0,sigp;
             0,0,0,0;
             0,0,0,0];

  
nti      = 20; 
nimp     = mN-1;
nt       = nti+nimp;
impinp   = [1,0,0,0];
inputs   = [zeros(nimp-1,4);impinp;zeros(nti,4)];
load cgg/cgg.dat
nt       = length(cgg);
inputs   = zeros(nt,4);
inputs(:,1) = cgg(:,2)/400;
[Xt,Xtl] = dlsim(G, H, G, H, inputs, [Xbar;1]);
pt       = exp(Xt(:,1));
wt       = exp(Xt(:,2));
Kit      = exp(Xt(:,3:N+2));
Kitl     = exp(Xtl(:,3:N+2));
Cjt      = exp(Xt(:,N+3:N+Nw+2)); 
Cjtl     = exp(Xtl(:,N+3:N+Nw+2)); 
Mdjt     = exp(Xt(:,N+Nw+3:N+2*Nw+2)); 
Yt       = exp(Xt(:,N+2*Nw+3));
rt       = Xt(:,N+2*Nw+4);
epsrt    = exp(Xtl(:,nZ+1));
Gt       = exp(Xtl(:,nZ+3));
At       = exp(Xtl(:,nZ+5));
pshkt    = exp(Xtl(:,nZ+7));

pvec     = [ps*ones(1,N-1)./mu.^[N-1:-1:1],pt']';
wvec     = [ws*ones(1,Nw-1)./mu.^[Nw-1:-1:1],wt']';
Pit      = zeros(nt,N);
Wit      = zeros(nt,Nw);
for i=1:N;
  Pit(:,i)  = pvec(N-i+1:nt-i+N);
end;
for i=1:Nw;
  Wit(:,i)  = wvec(Nw-i+1:nt-i+Nw);
end; 
Pbt      = (mean((Pit').^(-theta/thet1))').^(-thet1/theta);    
Wbt      = (mean((Wit').^(-v/v1))').^(-v1/v);
Pt       = Pbt.* mu.^[0:nt-1]';
Wt       = Wbt.* mu.^[0:nt-1]';
Yit      = ((Pt.^(1/thet1).*Yt)*ones(1,N))./Pit.^(1/thet1);
Ldit     = (Yit./(Kitl(:,[N,1:N-1]).^alpha)).^(1/alph1);
Ldit     = Ldit./(At*ones(1,N));
Ldt      = mean(Ldit')';            
KK       = Kit./Kitl(:,[N,1:N-1]);
Iit      = (1+b*delta-sqrt(1+2*b-2*b*KK))/b .* Kitl(:,[N,1:N-1]);     
It       = mean(Iit')';          
Lsjt     = ( (Ldt.*Wbt.^(1/v1))*ones(1,Nw) )./ Wit.^(1/v1);
Lt       = mean(Lsjt')';

%
% Plot variables of interest
%
t        = cgg(:,1);
dy       = cgg(:,3);
dpi      = cgg(:,4);
figure(1)
plot(t,100*log(Yt(1:nt)/ys),t,dy)
title('GDP Relative to Trend')
ylabel('Percent deviations')
legend('Model','Data')
figure(2)
plot(t,100*(([1;Pbt(2:nt)./Pbt(1:nt-1)]).^4-1),t,dpi)
title('Inflation Relative to Mean')
ylabel('Percent deviations')
legend('Model','Data')

