function R = res_niscmemsw_ms(Z,param)
%RES_NISCMEMSW_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-24-00
%                  Revised, ERM, 3-26-00 

%______________________________________________________________________________
%
% PARAMETERS


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);
v          = param(46);

alph1      = 1-alpha;
delt1      = 1-delta;
ome2       = 1-ome1;
rho1       = 1-rho;
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         = 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)'];
pf         = [X(2*N-2:-1:N);z(2,mN:-1:1)'];
pfa        = [X(3*N-3:-1:2*N-1);z(3,mN:-1:1)'];
pha        = [X(4*N-4:-1:3*N-2);z(4,mN:-1:1)'];
           % 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)'];
kai        = [X(5*N-3:6*N-4)';z(N+5:2*N+4,N:-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}]

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

w          = [X(nX-Nw+1:-1:6*N-3);z(2*N+7,mN:-1:1)'];
wa         = [X(nX:-1:6*N+Nw-4);z(2*N+8,mN:-1:1)'];
           % 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);   
mua        = Z(nz*mN+nX+4*mN-2:-1:nz*mN+nX+2*mN);
           % 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);
ga         = Z(nz*mN+nX+6*mN-2:-1:nz*mN+nX+5*mN-1);
           % 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)

kk         = ki(2:mN+1,:) ./ki(1:mN,[N,1:N-1]);
kka        = kai(2:mN+1,:)./kai(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]);
xai        = (1+b*delta-sqrt(1+2*b-2*b*kka))/b.* kai(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}]

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,:)  = [1;cumprod(mu(mN+i-2:-1:i+mN-N))]';
  pmuai(i,:) = [1;cumprod(mua(mN+i-2:-1:i+mN-N))]';
  wmui(i,:)  = [1;cumprod(mu(mN+i-2:-1:i+mN-Nw))]';
  wmuai(i,:) = [1;cumprod(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        = (sum( (phi' ./pmui') .^(-theta/thet1) )'/N).^(-thet1/theta);
pfb        = (sum( (pfi' ./pmui') .^(-theta/thet1) )'/N).^(-thet1/theta);
pfab       = (sum( (pfai'./pmuai').^(-theta/thet1) )'/N).^(-thet1/theta);
phab       = (sum( (phai'./pmuai').^(-theta/thet1) )'/N).^(-thet1/theta);
           % phb = [phb_{t-1}, ..., phb_{t+mN-2}]'

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

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

yhi        = (y .*(ome1*p) .^(1/rho1).*phb .^((theta-rho)/rho1/thet1)* ...
              o') .* (pmui ./phi).^(1/thet1);
yfi        = (y .*(ome2*p) .^(1/rho1).*pfb .^((theta-rho)/rho1/thet1)* ...
              o') .* (pmui ./pfi).^(1/thet1);
yfai       = (ya.*(ome1*pa).^(1/rho1).*pfab.^((theta-rho)/rho1/thet1)* ...
              o') .* (pmuai./pfai).^(1/thet1);
yhai       = (ya.*(ome2*pa).^(1/rho1).*phab.^((theta-rho)/rho1/thet1)* ...
              o') .* (pmuai./phai).^(1/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+N-1} yh_{2,t+N-1} ... yh_{N,t+N-1}]

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

Lsj        = (Ld .*wb .^(1/v1) * oj') .*(wmui ./wi) .^(1/v1);
Lsaj       = (Lda.*wab.^(1/v1) * oj') .*(wmuai./wai).^(1/v1);


mdj        = mu(mN:2*mN-1)./p     * oj';
mdaj       = mua(mN:2*mN-1)./pa   * oj';
cj         = (y-mean(xi')'-g)    * oj';
caj        = (ya-mean(xai')'-ga) * oj' ;

z0         = [cj(:);caj(:);mdj(:);mdaj(:)];
del        = diag(z0)*.0001;
for i=1:10;
  cj(:)        = z0(1:mN*Nw);
  caj(:)       = z0(mN*Nw+1:2*mN*Nw);
  mdj(:)       = z0(2*mN*Nw+1:3*mN*Nw);
  mdaj(:)      = z0(3*mN*Nw+1:4*mN*Nw);
 
  [U,UC,UM]    = util(cj(:),mdj(:),Lsj(:),param);
  [Ua,UCa,UMa] = util(caj(:),mdaj(:),Lsaj(:),param);

  res          = [mean(cj')'+mean(xi')'+g-y;
                  UC(mN+1:mN*Nw)-kron(ones(Nw-1,1),UC(1:mN));
                  mean(caj')'+mean(xai')'+ga-ya;
                  UCa(mN+1:mN*Nw)-kron(ones(Nw-1,1),UCa(1:mN));
                  mean(mdj')'-mu(mN:2*mN-1)./p;
                  UM(mN+1:mN*Nw)-kron(ones(Nw-1,1),UM(1:mN));
                  mean(mdaj')'-mua(mN:2*mN-1)./pa;
                  UMa(mN+1:mN*Nw)-kron(ones(Nw-1,1),UMa(1:mN))];

  for j=1:4*mN*Nw;
    zp           = z0 + del(:,j);
    cj(:)        = zp(1:mN*Nw);
    caj(:)       = zp(mN*Nw+1:2*mN*Nw);
    mdj(:)       = zp(2*mN*Nw+1:3*mN*Nw);
    mdaj(:)      = zp(3*mN*Nw+1:4*mN*Nw);

    [U,UC,UM]    = util(cj(:),mdj(:),Lsj(:),param);
    [Ua,UCa,UMa] = util(caj(:),mdaj(:),Lsaj(:),param);
    dres(:,j) ....
               = ([mean(cj')'+mean(xi')'+g-y;
                   UC(mN+1:mN*Nw)-kron(ones(Nw-1,1),UC(1:mN));
                   mean(caj')'+mean(xai')'+ga-ya;
                   UCa(mN+1:mN*Nw)-kron(ones(Nw-1,1),UCa(1:mN));
                   mean(mdj')'-mu(mN)/p(1);
                   UM(mN+1:mN*Nw)-kron(ones(Nw-1,1),UM(1:mN));
                   mean(mdaj')'-mua(mN)/pa(1);
                   UMa(mN+1:mN*Nw)-kron(ones(Nw-1,1),UMa(1:mN))]-res)/del(j,j);
  end;
  z0           = z0 - dres\res;

end;

cj(:)      = z0(1:mN*Nw);
caj(:)     = z0(mN*Nw+1:2*mN*Nw);
mdj(:)     = z0(2*mN*Nw+1:3*mN*Nw);
mdaj(:)    = z0(3*mN*Nw+1:4*mN*Nw);

[U,UC,UM,UL]     = util(cj,mdj,Lsj,param);
[Ua,UCa,UMa,ULa] = util(caj,mdaj,Lsaj,param);


rw         = wb./p;
rwa        = wab./pa;
           % 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);
mca        = rwa(1:N)./(alph1*(diag(kai(1:N,[N,1:N-1]))./ ...
             diag(Ldai(1:N,:))).^alpha.*aa(1:N).^alph1);
           % mc = [mc_{1,t}, mc_{2,t+1}, ...mc_{N,t+N-1}]'

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

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

bm1        = beta.^[0:N-1]' .* ([1;cumprod(mu(mN:mN+N-2))]).^(1/thet1);
bm2        = beta.^[0:N-1]' .* ([1;cumprod(mu(mN:mN+N-2))]).^(theta/thet1);
bma1       = beta.^[0:N-1]' .* ([1;cumprod(mua(mN:mN+N-2))]).^(1/thet1);
bma2       = beta.^[0:N-1]' .* ([1;cumprod(mua(mN:mN+N-2))]).^(theta/thet1);
uc         = UC(1:N,1);
uca        = UCa(1:N,1);
 
R(1,1)     = ph(N)*theta *sum(uc ./p(1:N).*lh .*bm2)  ...
                         -sum(uc .*mc .*lh .*bm1);
R(2,1)     = pf(N)*theta *sum(uca./(p(1:N).*q).*lf .*bm2) ...
                         -sum(uca.*mca.*lf .*bm1);
R(3,1)     = pfa(N)*theta*sum(uca./pa(1:N)    .*lfa.*bma2) ...
                         -sum(uca.*mca.*lfa.*bma1);
R(4,1)     = pha(N)*theta*sum(uc .*q./pa(1:N) .*lha.*bma2) ...
                         -sum(uc .*mc .*lha.*bma1);

%
% 2. 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(5:N+4,1) = UC(1,1)*o./(1-dphi) - beta*UC(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));

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

R(N+5:2*N+4,1)  ...
           = UCa(1)*o./(1-dphi) - beta*UCa(2)*(rwa(2)*alpha/alph1* ...
                Ldai(2,[2:N,1])'./kai(2,:)'+(1-delta- ...
                phin+dphin.*xai(2,[2:N,1])'./kai(2,:)')./(1-dphin));

%
% 3. MONEY DEMAND EQUATIONS
%

R(2*N+5,1) = UM(1,1)/UC(1,1)  +beta*UC(2,2)/UC(1,1)*p(1)/p(2)/mu(mN)-1;
R(2*N+6,1) = UMa(1,1)/UCa(1,1)+beta*UCa(2,2)/UCa(1,1)*pa(1)/pa(2)/mua(mN)-1;


%
% 4. WAGE EQUATION
%

bm1        = beta.^[0:Nw-1]' .* ([1;cumprod(mu(mN:mN+Nw-2))]).^(1/v1);
bm2        = beta.^[0:Nw-1]' .* ([1;cumprod(mu(mN:mN+Nw-2))]).^(v/v1);
bma1       = beta.^[0:Nw-1]' .* ([1;cumprod(mua(mN:mN+Nw-2))]).^(1/v1);
bma2       = beta.^[0:Nw-1]' .* ([1;cumprod(mua(mN:mN+Nw-2))]).^(v/v1);
wb         = wb(1:Nw);
wab        = wab(1:Nw);
p          = p(1:Nw);
pa         = pa(1:Nw);
uc         = diag(UC(1:Nw,:));
uca        = diag(UCa(1:Nw,:));
ul         = diag(UL(1:Nw,:));
ula        = diag(ULa(1:Nw,:));

R(2*N+7,1) = w(Nw)*v *sum(uc ./p .*wb .^(1/v1).*Ld .*bm2) ...
                     +sum(ul     .*wb .^(1/v1).*Ld .*bm1);
R(2*N+8,1) = wa(Nw)*v*sum(uca./pa.*wab.^(1/v1).*Lda.*bma2) ...
                     +sum(ula    .*wab.^(1/v1).*Lda.*bma1);


