function R = res_iscmemsw_ms(Z,param)
%RES_ISCMEMSW_MS Residuals for the international sticky price model with 
%                complete markets, exogenous growth rates for money, 
%                sticky wages, multiple shocks, and N>1.

%                Ellen R. McGrattan, 3-23-00
%                Revised, ERM, 4-27-00

%______________________________________________________________________________
%
% PARAMETERS


ns         = 50;
N          = param(2);
Nw         = param(3);
beta       = param(11);
ome1       = param(12);
rho        = param(13);
theta      = param(14);
alpha      = param(15);
delta      = param(16);
b          = param(17);
mus        = param(19);
gs         = param(36);
v          = param(46);
phs        = param(ns);
ys         = param(ns+1);
ws         = param(ns+2);
cjs        = param(ns+3:ns+2+Nw);
mdjs       = param(ns+3+Nw:ns+2+2*Nw);
Lsjs       = param(ns+3+2*Nw:ns+2+3*Nw);
UC         = param(ns+3+3*Nw:ns+2+4*Nw);
UM         = param(ns+3+4*Nw:ns+2+5*Nw);
UL         = param(ns+3+5*Nw:ns+2+6*Nw);
UCC        = param(ns+3+6*Nw:ns+2+7*Nw);
UMM        = param(ns+3+7*Nw:ns+2+8*Nw);
ULL        = param(ns+3+8*Nw:ns+2+9*Nw);
UCM        = param(ns+3+9*Nw:ns+2+10*Nw);
UCL        = param(ns+3+10*Nw:ns+2+11*Nw);
ULM        = param(ns+3+11*Nw:ns+2+12*Nw);
kis        = param(ns+3+12*Nw:ns+2+12*Nw+N);
xis        = param(ns+3+12*Nw+N:ns+2+12*Nw+2*N);
yhis       = param(ns+3+12*Nw+2*N:ns+2+12*Nw+3*N);
yhais      = param(ns+3+12*Nw+3*N:ns+2+12*Nw+4*N);
Ldis       = param(ns+3+12*Nw+4*N:ns+2+12*Nw+5*N);
mcis       = param(ns+3+12*Nw+5*N:ns+2+12*Nw+6*N);
alph1      = 1-alpha;
delt1      = 1-delta;
ome2       = 1-ome1;
rho1       = 1-rho;
thet1      = 1-theta;
v1         = 1-v;
phis       = b/2*(xis./kis([N,1:N-1])-delta).^2;
dphis      = b*(xis./kis([N,1:N-1])-delta);
d2phis     = b*ones(N,1);
Lds        = sum(Ldis)/N;
UC         = UC(1);
UM         = UM(1);
ps         = mus/mean(mdjs);
mvec1      = mus.^(theta/thet1*[0:N-1]);
mvec2      = mus.^(v/v1*[0:Nw-1]);
phbs       = ps/(ome1^(1/rho1)+ome2^(1/rho1))^(-rho1/rho);
rws        = ws/ps * mean((mus^(v/v1)).^[0:Nw-1]')^(-v1/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
%

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

ph         = [X(N-1:-1:1);z(1,mN:-1:1)'] - log(phs);
pf         = [X(2*N-2:-1:N);z(2,mN:-1:1)'] - log(phs);
pfa        = [X(3*N-3:-1:2*N-1);z(3,mN:-1:1)'] - log(phs);
pha        = [X(4*N-4:-1:3*N-2);z(4,mN:-1:1)'] - log(phs);
           % ph  = [ph_{t-N}, ...,ph_t, ... ph_{t+mN-2}]'   (N+mN-1 x 1)
           % pf  = [pf_{t-N}, ...,pf_t, ... pf_{t+mN-2}]'   (N+mN-1 x 1)
           % pf* = [pf*_{t-N},...,pf*_t, ...pf*_{t+mN-2}]'  (N+mN-1 x 1)
           % ph* = [ph*_{t-N},...,ph*_t, ...ph*_{t+mN-2}]'  (N+mN-1 x 1)

ki         = [X(4*N-3:5*N-4)';z(5:N+4,mN:-1:1)']    - log(o1*kis');
kai        = [X(5*N-3:6*N-4)';z(N+5:2*N+4,mN:-1:1)'] - log(o1*kis');
           % 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}]

y          = z(2*N+5,mN:-1:1)' - log(ys);
ya         = z(2*N+6,mN:-1:1)' - log(ys);
           % y   = [y_t, ... y_{t+mN-1}]'

w          = [X(nX-Nw+1:-1:6*N-3);z(2*N+7,mN:-1:1)'] - log(ws);
wa         = [X(nX:-1:6*N+Nw-4);z(2*N+8,mN:-1:1)']   - log(ws);
           % w   = [w_{t-Nw}, ...,w_t, ... w_{t+mN-1}]'  (Nw+mN-1 x 1)

mu         = Z(nz*mN+nX+2*mN-1:-1:nz*mN+nX+1)    - log(mus);
mua        = Z(nz*mN+nX+4*mN-2:-1:nz*mN+nX+2*mN) - log(mus);
           % mu  = [mu_{t-mN+1}, ..., mu_{t+mN-1}]'   (2*mN-1 x 1)

g          = Z(nz*mN+nX+5*mN-2:-1:nz*mN+nX+4*mN-1) - log(gs);
ga         = Z(nz*mN+nX+6*mN-2:-1:nz*mN+nX+5*mN-1) - log(gs);
           % g   = [g_{t-mN+1},  ..., g_{t+mN-1}]'   (mN x 1)


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

tem        = o*(kis./(xis.*(1-dphis)))';
xi         = ki(1:mN,[N,1:N-1])+(ki(2:mN+1,:)-ki(1:mN,[N,1:N-1])) .* tem; 
xai        = kai(1:mN,[N,1:N-1])+(kai(2:mN+1,:)-kai(1:mN,[N,1:N-1])) .* tem; 
           % 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}]


phi        = zeros(mN,N);
pfi        = zeros(mN,N);
pfai       = zeros(mN,N);
phai       = zeros(mN,N);
pmui       = zeros(mN,N);
pmuai      = zeros(mN,N);
wmui       = zeros(mN,Nw);
wmuai      = zeros(mN,Nw);
for i=1:mN;
  phi(i,:)   = ph(i+N-1:-1:i)'; 
  pfi(i,:)   = pf(i+N-1:-1:i)'; 
  pfai(i,:)  = pfa(i+N-1:-1:i)'; 
  phai(i,:)  = pha(i+N-1:-1:i)'; 
  wi(i,:)    = w(i+Nw-1:-1:i)'; 
  wai(i,:)   = wa(i+Nw-1:-1:i)'; 
  pmui(i,:)  = [0;cumsum(mu(mN+i-2:-1:i+mN-N))]';
  pmuai(i,:) = [0;cumsum(mua(mN+i-2:-1:i+mN-N))]';
  wmui(i,:)  = [0;cumsum(mu(mN+i-2:-1:i+mN-Nw))]';
  wmuai(i,:) = [0;cumsum(mua(mN+i-2:-1:i+mN-Nw))]';
end;
           % phi = [ph_{t-1},   ph_{t-2},    ..., ph_{t-N},
           %        ph_t    ,   ph_{t-1},    ..., ph_{t-N+1},
           %           :            :             :
           %        ph_{t+mN-2},ph_{t+mN-3}, ..., ph_{t+mN-N-1}]
           %
           % pmui= [0, mu_{t-1},   ...,mu_{t-1}+mu_{t-2}+...+mu_{t-N+1},
           %        0, mu_t,       ...,mu_t    +mu_{t-1}+...+mu_{t-N+2},
           %           :               :             
           %        0, mu_{t+mN-2},...,mu_{t+mN-2}+mu_{t+mN-3}...+mu_{t+mN-N}]
           %
           % wi  = [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}]
           %
           % wmui= [0, mu_{t-1},   ...,mu_{t-1}+mu_{t-2}+...+mu_{t-Nw+1},
           %        0, mu_t,       ...,mu_t    +mu_{t-1}+...+mu_{t-Nw+2},
           %           :               :             
           %        0, mu_{t+mN-2},...,mu_{t+mN-2}+mu_{t+mN-3}...+mu_{t+mN-Nw}]
           %


phb        = (phi-pmui)*mvec1'/sum(mvec1);
pfb        = (pfi-pmui)*mvec1'/sum(mvec1);
pfab       = (pfai-pmuai)*mvec1'/sum(mvec1);
phab       = (phai-pmuai)*mvec1'/sum(mvec1);
           % phb = [phb_{t-1},phb_t, ..., phb_{t+mN-2}]'

wb         = (wi-wmui)*mvec2'/sum(mvec2);
wab        = (wai-wmuai)*mvec2'/sum(mvec2);
           % wb  = [wb_t, wb_{t+1},..., wb_{t+mN-1}]'

p          = (ps/phbs)^(rho/rho1)*(ome1^(1/rho1)*phb+ome2^(1/rho1)*pfb);
pa         = (ps/phbs)^(rho/rho1)*(ome1^(1/rho1)*pfab+ome2^(1/rho1)*phab);
           % p   = [p_t,p_{t+1}, ..., p_{t+mN-1}]'

rw         = wb-p;
rwa        = wab-pa;
           % rw = [rw_t,rw_{t+1},...rw_{t+mN-1}]'

yhi        = (y+p/rho1-(rho-theta)*phb/rho1/thet1)*oi'-(phi-pmui)/thet1;
yfi        = (y+p/rho1-(rho-theta)*pfb/rho1/thet1)*oi'-(pfi-pmui)/thet1;
yfai       = (ya+pa/rho1-(rho-theta)*pfab/rho1/thet1)*oi'-(pfai-pmuai)/thet1;
yhai       = (ya+pa/rho1-(rho-theta)*phab/rho1/thet1)*oi'-(phai-pmuai)/thet1;
           % yhi = [yh_{1,t}      yh_{2,t}      ... yh_{N,t}
           %        yh_{1,t+1}    yh_{2,t+1}    ... yh_{N,t+1}
           %           :             :                 :
           %        yh_{1,t+mN-1} yh_{2,t+mN-1} ... yh_{N,t+mN-1}]

Ldi        = (yhi.*(o*yhis') + yhai.*(o*yhais'))./(alph1*o*(yhis'+yhais'))- ...
             alpha*ki(1:mN,[N,1:N-1])/alph1 - a*ones(1,N);
Ldai       = (yfai.*(o*yhis')+ yfi.*(o*yhais'))./(alph1*o*(yhis'+yhais'))- ...
             alpha*kai(1:mN,[N,1:N-1])/alph1 - aa*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         = Ldi* Ldis/sum(Ldis);
Lda        = Ldai*Ldis/sum(Ldis);
           % Ld  = [Ld_t,...,Ld_{t+mN-1}]' 

Lsj        = (Ld+wb/v1)*oj'+(wmui-wi)/v1;
Lsaj       = (Lda+wab/v1)*oj'+(wmuai-wai)/v1;
           % Lsj  = [Ls_{1,t}       Ls_{2,t}          ... Ls_{Nw,t}
           %         Ls_{1,t+1}     Ls_{2,t+1}        ... Ls_{Nw,t+1}
           %            :              :                     :
           %         Ls_{1,t+mN-1}  Ls_{2,t+mN-1}     ... Ls_{Nw,t+mN-1}]


A          = [cjs'/Nw,zeros(1,Nw);
              -UCC(1)*cjs(1)*ones(Nw-1,1),diag(UCC(2:Nw).*cjs(2:Nw)), ...
              -UCM(1)*mdjs(1)*ones(Nw-1,1),diag(UCM(2:Nw).*mdjs(2:Nw));
              zeros(1,Nw),ones(1,Nw)/Nw;
              -UCM(1)*cjs(1)*ones(Nw-1,1),diag(UCM(2:Nw).*cjs(2:Nw)), ...
              -UMM(1)*mdjs(1)*ones(Nw-1,1),diag(UMM(2:Nw).*mdjs(2:Nw))];

for i=1:mN;
  tem      = A\[ys*y(i)-gs*g(i)-xi(i,:)*xis/N;
                UCL(1)*Lsjs(1)*Lsj(i,1)-UCL(2:Nw).*Lsjs(2:Nw).*Lsj(i,2:Nw)';
                mu(i+mN-1)-p(i);
                ULM(1)*Lsjs(1)*Lsj(i,1)-ULM(2:Nw).*Lsjs(2:Nw).*Lsj(i,2:Nw)'];
  cj(i,:)  = tem(1:Nw)';
  mdj(i,:) = tem(Nw+1:2*Nw)';
  tem      = A\[ys*ya(i)-gs*ga(i)-xai(i,:)*xis/N;
                UCL(1)*Lsjs(1)*Lsaj(i,1)-UCL(2:Nw).*Lsjs(2:Nw).*Lsaj(i,2:Nw)';
                mua(i+mN-1)-pa(i);
                ULM(1)*Lsjs(1)*Lsaj(i,1)-ULM(2:Nw).*Lsjs(2:Nw).*Lsaj(i,2:Nw)'];
  caj(i,:) = tem(1:Nw)';
  mdaj(i,:)= tem(Nw+1:2*Nw)';
end;
           % cj   = [c_{1,t}        c_{2,t}           ... c_{Nw,t}
           %         c_{1,t+1}      c_{2,t+1}         ... c_{Nw,t+1}
           %           :              :                     :
           %         c_{1,t+mN-1}   c_{2,t+mN-1}      ... c_{Nw,t+mN-1}]
           %
           % mdj  = [md_{1,t}       md_{2,t}          ... md_{Nw,t}
           %         md_{1,t+1}     md_{2,t+1}        ... md_{Nw,t+1}
           %            :              :                     :
           %         md_{1,t+mN-1}  md_{2,t+mN-1}     ... md_{Nw,t+mN-1}]

ucj        = cj  .* ( o* (UCC.*cjs)' ) + Lsj .* ( o* (UCL.*Lsjs)' ) + ...
             mdj .* ( o* (UCM.*mdjs)'); 
ucaj       = caj .* ( o* (UCC.*cjs)' ) + Lsaj.* ( o* (UCL.*Lsjs)' ) + ...
             mdaj.* ( o* (UCM.*mdjs)'); 
ulj        = cj  .* ( o* (UCL.*cjs)' ) + Lsj .* ( o* (ULL.*Lsjs)' ) + ...
             mdj .* ( o* (ULM.*mdjs)'); 
ulaj       = caj .* ( o* (UCL.*cjs)' ) + Lsaj.* ( o* (ULL.*Lsjs)' ) + ...
             mdaj.* ( o* (ULM.*mdjs)'); 
           % ucj  = [Uc_{1,t}       Uc_{2,t}          ... Uc_{Nw,t}
           %         Uc_{1,t+1}     Uc_{2,t+1}        ... Uc_{Nw,t+1}
           %            :              :                     :
           %         Uc_{1,t+mN-1}  Uc_{2,t+mN-1}     ... Uc_{Nw,t+1}

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

q          = (ucaj(1:N,1)-ucj(1:N,1))/UC;
           % q  = [q_t,q_{t+1}, ...q_{t+N-1}]'

lh         = p(1:N)/rho1+phb(1:N)*(theta-rho)/rho1/thet1+y(1:N);
lf         = p(1:N)/rho1+pfb(1:N)*(theta-rho)/rho1/thet1+y(1:N);
lfa        = pa(1:N)/rho1+pfab(1:N)*(theta-rho)/rho1/thet1+ya(1:N);
lha        = pa(1:N)/rho1+phab(1:N)*(theta-rho)/rho1/thet1+ya(1:N);
           % lh = [lam_{h,t},lam_{h,t+1}, ...lam_{h,t+N-1}]'

uc         = ucj(1:N,1);
uca        = ucaj(1:N,1);

%______________________________________________________________________________
%
% LINEARIZED RESIDUALS:
%
% 1. PRICING EQUATIONS
%

w1             = (beta*mus^(1/thet1)).^[0:N-1]'.*mcis; w1 = w1/sum(w1);
w2             = (beta*mus^(theta/thet1)).^[0:N-1]';   w2 = w2/sum(w2);
pmui           = diag(pmui(1:N,:));
pmuai          = diag(pmuai(1:N,:));

R(1,1)         = ph(N) -w1'*(uc/UC  +lh  +mc       +pmui/thet1) ...
                       +w2'*(uc/UC  +lh  -p(1:N)   +pmui*theta/thet1);
R(2,1)         = pf(N) -w1'*(uca/UC +lf  +mca      +pmui/thet1) ...
                       +w2'*(uca/UC +lf  -q-p(1:N) +pmui*theta/thet1);
R(3,1)         = pfa(N)-w1'*(uca/UC +lfa +mca      +pmuai/thet1) ...
                       +w2'*(uca/UC +lfa -pa(1:N)  +pmuai*theta/thet1);
R(4,1)         = pha(N)-w1'*(uc/UC  +lha +mc       +pmuai/thet1) ...
                       +w2'*(uc/UC  +lha +q-pa(1:N)+pmuai*theta/thet1);

%
% 2. MONOPOLISTS' EULER EQUATIONS FOR CAPITAL
%

tem1           = beta*rws*alpha/alph1*(1-dphis).*Ldis([2:N,1])./kis;
tem2           = d2phis.*xis./kis([N,1:N-1])./(1-dphis);
tem3           = d2phis([2:N,1]).*xis([2:N,1])./kis./(1-dphis([2:N,1]));
tem4           = (1-tem1) .* tem3 + ...
                 beta*(1-dphis).*d2phis([2:N,1]).*  ...
                 (xis([2:N,1])./kis).^2./(1-dphis([2:N,1]));

R(5:N+4,1)     =-((ucj(2,2)-ucj(1,1))/UC ...
                 +tem1 .* (rw(2)+Ldi(2,[2:N,1])'-ki(2,:)') ...
                 -tem2 .* (xi(1,:)'-ki(1,[N,1:N-1])') ...
                 +tem4 .* (xi(2,[2:N,1])'-ki(2,:)'));

R(N+5:2*N+4,1) =-((ucaj(2,2)-ucaj(1,1))/UC ...
                 +tem1 .* (rwa(2)+Ldai(2,[2:N,1])'-kai(2,:)') ...
                 -tem2 .* (xai(1,:)'-kai(1,[N,1:N-1])') ...
                 +tem4 .* (xai(2,[2:N,1])'-kai(2,:)'));

%
% 3. MONEY DEMAND EQUATIONS
%

R(2*N+5,1)     = UM/UC*((UCM(1)/UM-UCC(1)/UC)*cjs(1) *cj(1,1)+ ...
                        (ULM(1)/UM-UCL(1)/UC)*Lsjs(1)*Lsj(1,1)+ ...
                        (UMM(1)/UM-UCM(1)/UC)*mdjs(1)*mdj(1,1))+ ...
               beta/mus*((ucj(2,2)-ucj(1,1))/UC+ ...
                         p(1)-p(2)-mu(mN));

R(2*N+6,1)     = UM/UC*((UCM(1)/UM-UCC(1)/UC)*cjs(1) *caj(1,1)+ ...
                        (ULM(1)/UM-UCL(1)/UC)*Lsjs(1)*Lsaj(1,1)+ ...
                        (UMM(1)/UM-UCM(1)/UC)*mdjs(1)*mdaj(1,1))+ ...
               beta/mus*((ucaj(2,2)-ucaj(1,1))/UC+ ...
                         pa(1)-pa(2)-mua(mN));

%
% 4. WAGE EQUATION
%

w1             = (beta*mus^(1/v1)).^[0:Nw-1]'.*UL; w1 = w1/sum(w1);
w2             = (beta*mus^(v/v1)).^[0:Nw-1]';     w2 = w2/sum(w2);
wmui           = diag(wmui(1:Nw,:));  
wmuai          = diag(wmuai(1:Nw,:));
ulj            = diag(ulj(1:Nw,:));
ucj            = diag(ucj(1:Nw,:));
ulaj           = diag(ulaj(1:Nw,:));
ucaj           = diag(ucaj(1:Nw,:));
Ld             = Ld(1:Nw);
Lda            = Lda(1:Nw);
wb             = wb(1:Nw);
wab            = wab(1:Nw);
p              = p(1:Nw);
pa             = pa(1:Nw);

R(2*N+7,1)     = w(Nw) -w1'*(ulj ./UL+wb/v1 +Ld +wmui/v1) ...
                       +w2'*(ucj  /UC+wb/v1 +Ld +wmui*v/v1-p);
R(2*N+8,1)     = wa(Nw)-w1'*(ulaj./UL+wab/v1+Lda+wmuai/v1) ...
                       +w2'*(ucaj /UC+wab/v1+Lda+wmuai*v/v1-pa);
