function R = res_triple(Z,param)
%RES_GALI    Residuals for a sticky price/sticky wage model with 
%            staggered price and wage setting, complete markets, 
%            habit persistence in preferences, a Taylor rule for 
%            the Fed, and multiple shocks.

%            Ellen R. McGrattan, 3-24-00
%            Revised, ERM, 3-20-04

%______________________________________________________________________________
%
% PARAMETERS

N          = param(2);
Nw         = param(3);
beta       = param(11);
alpha      = param(12);
delta      = param(13);
b          = param(14);
theta      = param(15);
v          = param(16);
mu         = param(17);
irule      = param(19:28);
rs         = mu/beta-1;
ys         = param(31+N+2*Nw);

alph1      = 1-alpha;
delt1      = 1-delta;
thet1      = 1-theta;
v1         = 1-v;

mN         = max(N,Nw);
oi         = ones(N,1);
oj         = ones(Nw,1);
o          = ones(mN,1);
o1         = ones(mN+1,1);

%______________________________________________________________________________
%
% VARIABLES APPEARING IN THE RESIDUALS
%

Z          = exp(Z);
nz         = N+2*Nw+4;
nX         = 2*N+2*Nw+6;
z          = zeros(nz,mN);
z(:)       = Z(1:nz*mN);
z(nz,:)    = log(z(nz,:));
X          = Z(nz*mN+1:nz*mN+nX);

p          = [X(N-1:-1:1);z(1,mN:-1:1)'];
           % p  = [p_{t-N}, ...,p_t, ... p_{t+mN-2}]'   (N+mN-1 x 1)

w          = [X(N+Nw+1:-1:N+3);z(2,mN:-1:1)'];
           % w   = [w_{t-Nw}, ...,w_t, ... w_{t+mN-1}]'  (Nw+mN-1 x 1)

ki         = [X(N+Nw+2:2*N+Nw+1)';z(3:N+2,mN:-1:1)'];
           % ki  = [k_{1,t-1},     k_{2,t-1},   ...   k_{N,t-1}
           %        k_{1,t}        k_{2,t},     ...   k_{N,t}
           %           :             :                 :
           %        k_{1,t+mN-1},  k_{2,t+mN-1}, ...  k_{N,t+mN-1}]

cj         = [X(2*N+Nw+2:2*N+2*Nw+1)';z(N+3:N+Nw+2,mN:-1:1)'];
           % cj  = [c_{1,t-1},     c_{2,t-1},    ...  c_{Nw,t-1}
           %        c_{1,t}        c_{2,t},      ...  c_{Nw,t}
           %           :             :                 :
           %        c_{1,t+mN-1},  c_{2,t+mN-1}, ...  c_{Nw,t+mN-1}]

mdj        = z(N+Nw+3:N+2*Nw+2,mN:-1:1)';
           % mdj  =[md_{1,t}       md_{2,t},      ... md_{Nw,t}
           %           :             :                :
           %        md_{1,t+mN-1}, md_{2,t+mN-1}, ... md_{Nw,t+mN-1}]


y          = z(N+2*Nw+3,mN:-1:1)';
yl1        = X(2*N+2*Nw+2);
yl2        = X(2*N+2*Nw+3);
           % y   = [y_t, ... y_{t+mN-1}]'
           % yli = y_{t-i}

r          = z(N+2*Nw+4,mN:-1:1)';
rl1        = log(X(2*N+2*Nw+4));
rl2        = log(X(2*N+2*Nw+5));
rl3        = log(X(2*N+2*Nw+6));
           % r   = [r_t, ... r_{t+mN-1}]'
           % rli = r_{t-i}

eps        = log(Z(nz*mN+nX+1));
           % eps  = eps_t

g          = Z(nz*mN+nX+mN+1:-1:nz*mN+nX+2);
           % g   = [g_t,  ..., g_{t+mN-1}]'   (mN x 1)

a          = Z(nz*mN+nX+2*mN+1:-1:nz*mN+nX+mN+2);
           % a   = [a_t,  ..., a_{t+mN-1}]'   (mN x 1)

pshk       = Z(nz*mN+nX+3*mN+1:-1:nz*mN+nX+2*mN+2);
           % pshk= [pshk_t,  ..., pshk_{t+mN-1}]'   (mN x 1)

kk         = ki(2:mN+1,:) ./ki(1:mN,[N,1:N-1]);
xi         = (1+b*delta-sqrt(1+2*b-2*b*kk))/b .* ki(1:mN,[N,1:N-1]);
           % x   = [x_{1,t}       x_{2,t}       ... x_{N,t}
           %        x_{1,t+1}     x_{2,t+1}     ... x_{N,t+1}
           %           :             :                 :
           %        x_{1,t+mN-1}  x_{2,t+mN-1}  ... x_{N,t+mN-1}]

pi         = zeros(mN,N);
mu1        = mu.^[0:N-1]';
mui        = o*mu1';
muj        = o*(mu.^[0:Nw-1]')';

for i=1:mN;
  pi(i,:)    = p(i+N-1:-1:i)';
  wj(i,:)    = w(i+Nw-1:-1:i)';
end;
           % pi  = [p_{t-1},   p_{t-2},    ..., p_{t-N},
           %        p_t    ,   p_{t-1},    ..., p_{t-N+1},
           %           :            :             :
           %        p_{t+mN-2},p_{t+mN-3}, ..., p_{t+mN-N-1}]
           %
           % wj  = [w_{t-1},    w_{t-2},     ..., w_{t-Nw},
           %        w_t    ,    w_{t-1},     ..., w_{t-Nw+1},
           %           :            :             :
           %        w_{t+mN-2}, w_{t+mN-3},  ..., w_{t+mN-Nw-1}]
           %
pil1       = [pi(1,2:N),X(N)];
pil2       = [pil1(1,2:N),X(N+1)];
pil3       = [pil2(1,2:N),X(N+2)];
           % pili = [p_{t-i-1}, ..., p_{t-N-i}]

pb         = (sum( (pi' ./mui') .^(-theta/thet1) )'/N).^(-thet1/theta);
           % pb = [pb_{t-1}, ..., pb_{t+mN-2}]'

pbl1       = (sum( (pil1 ./mu1') .^(-theta/thet1) )'/N).^(-thet1/theta);
pbl2       = (sum( (pil2 ./mu1') .^(-theta/thet1) )'/N).^(-thet1/theta);
pbl3       = (sum( (pil3 ./mu1') .^(-theta/thet1) )'/N).^(-thet1/theta);
           % pbli = pb_{t-i-1}
                                                                                
wb         = (sum( (wj' ./muj') .^(-v/v1) )'/Nw).^(-v1/v);
           % wb  = [wb_t, wb_{t+1},..., wb_{t+mN-1}]'

yi         = (y .*  pb.^(1/thet1)* o') .* (mui ./pi).^(1/thet1);

           % yi = [y_{1,t}     y_{2,t}     ... y_{N,t}
           %       y_{1,t+1}   y_{2,t+1}   ... y_{N,t+1}
           %           :            :             :
           %       y_{1,t+N-1} y_{2,t+N-1} ... y_{N,t+N-1}]

Ldi        = ( yi./ki(1:mN,[N,1:N-1]).^alpha ).^(1/alph1) ./(a*ones(1,N));
           % Ldi =  [Ld_{1,t}       Ld_{2,t}       ... Ld_{N,t}
           %         Ld_{1,t+1}     Ld_{2,t+1}     ... Ld_{N,t+1}
           %            :              :                  :
           %         Ld_{1,t+mN-1}  Ld_{2,t+mN-1}  ... Ld_{N,t+mN-1}]

Ld         = sum(Ldi')'/N;
           % Ld  = [Ld_t,...,Ld_{t+mN-1}]'

Lsj        = (Ld .*wb .^(1/v1) * oj') .*(muj ./wj) .^(1/v1);

rw         = wb./pb;
           % rw = [rw_t,rw_{t+1},...rw_{t+mN-1}]'

mc         = rw(1:N)./(alph1*(diag(ki(1:N,[N,1:N-1]))./ ...
             diag(Ldi(1:N,:))).^alpha.*a(1:N).^alph1);
           % mc = [mc_{1,t}, mc_{2,t+1}, ...mc_{N,t+N-1}]'

lam        = pb.^(1/thet1).*y;
           % lam = [lam_t,lam_{t+1}, ...lam_{t+mN-1}]'

%[U,U1,U2,UM,UL]  = util(cj(2:mN+1,:),cj(1:mN,[Nw,1:Nw-1]),mdj,Lsj, ...
%                        pshk*ones(1,Nw),param);
[U,U1,U2,UM,UL] = util(cj(2:mN+1,:),mean(cj(1:mN,:)')'*ones(1,Nw), ...
                       mdj,Lsj,pshk*ones(1,Nw),param); %external

zetj       = UM .* ((1+r)./r * oj');
%______________________________________________________________________________
%
% LINEARIZED RESIDUALS:
%
% 1. PRICING EQUATION
%

bm1        = (beta*mu^(1/thet1)).^[0:N-1]';
bm2        = (beta*mu^(theta/thet1)).^[0:N-1]';
if (Nw>N);
  zet      = diag(zetj(1:N,1:N));
else
  j          = 1;
  for i=1:N;
    zet(i) = zetj(i,j);
    if j==Nw; j=1; end;
  end;
end;
zet = zet(:);
R(1,1)     = p(N)*theta *sum(zet ./pb(1:N) .*lam .*bm2)  ...
                        -sum(zet .*mc .*lam .*bm1);

%
% 2. WAGE EQUATION
%

bm1        = (beta*mu^(1/v1)).^[0:Nw-1]'; 
bm2        = (beta*mu^(v/v1)).^[0:Nw-1]';
zet        = diag(zetj(1:Nw,:));
ul         = diag(UL(1:Nw,:));

R(2,1)   = w(Nw)*v *sum(zet./pb(1:Nw) .*wb(1:Nw) .^(1/v1).*Ld(1:Nw).*bm2) ...
                   +sum(ul            .*wb(1:Nw) .^(1/v1).*Ld(1:Nw).*bm1);

%
% 3. MONOPOLISTS' EULER EQUATIONS FOR CAPITAL
%

phin       = b/2*(xi(2,[2:N,1])'./ki(2,:)'-delta).^2; 
dphi       = b*(xi(1,:)'./ki(1,[N,1:N-1])'-delta); 
dphin      = b*(xi(2,[2:N,1])'./ki(2,:)'-delta); 

R(3:N+2,1) = zetj(1,1)*o./(1-dphi) - beta*zetj(2,2)*(rw(2)*alpha/alph1* ...
                Ldi(2,[2:N,1])'./ki(2,:)'+(1-delta- ...
                phin+dphin.*xi(2,[2:N,1])'./ki(2,:)')./(1-dphin));

%
% 5. CONSUMPTION EQUATIONS (with internal/external habit)
%
%R(N+3:N+Nw+2,1)   = (U1(1,[Nw,1:Nw-1])+beta*U2(2,:)-zetj(1,[Nw,1:Nw-1]))';
R(N+3:N+Nw+2,1)  = (U1(1,[Nw,1:Nw-1])-zetj(1,[Nw,1:Nw-1]))'; %external

%
% 6. MONEY DEMAND EQUATIONS
%

R(N+Nw+3:N+2*Nw+2,1)   = (UM(1,[Nw,1:Nw-1])-zetj(1,[Nw,1:Nw-1])+ ....
                          beta*zetj(2,:)*pb(1)/(pb(2)*mu))';

%
% 7. RESOURCE CONSTRAINT
%
tem             = mean(cj(2:mN+1,:)')'+mean(xi')'+g-y;
R(N+2*Nw+3,1)   = tem(1);

%
% 8. INTEREST RATE EQUATION
%

R(N+2*Nw+4,1)   = r(1) -rs-irule' *[rl1-rs;rl2-rs;rl3-rs;
                                    log(pb(2))-log(pb(1));
                                    log(pb(1))-log(pbl1);
                                    log(pbl1)-log(pbl2);
                                    log(pbl2)-log(pbl3);
                                    log(y(1))-log(ys);
                                    log(yl1)-log(ys);
                                    log(yl2)-log(ys)]-eps;

