function R = res_iscmem(Z,param)
%RES_ISCMEM    Residuals for the international sticky price model with 
%              exogenous growth rates for money and N>1.

%              Ellen R. McGrattan, 10-11-96
%              Revised, ERM, 4-27-00

%______________________________________________________________________________
%
% PARAMETERS


ns         = 50;
N          = param(2);
beta       = param(11);
ome1       = param(12);
rho        = param(13);
theta      = param(14);
alpha      = param(15);
delta      = param(16);
b          = param(17);
mus        = param(19);
phs        = param(ns);
ys         = param(ns+1);
cs         = param(ns+2);
MPs        = param(ns+3);
UC         = param(ns+4);
UM         = param(ns+5);
UL         = param(ns+6);
UCC        = param(ns+7);
UMM        = param(ns+8);
ULL        = param(ns+9);
UCM        = param(ns+10);
UCL        = param(ns+11);
ULM        = param(ns+12);
kis        = param(ns+13:ns+12+N);
xis        = param(ns+13+N:ns+12+2*N);
yhis       = param(ns+13+2*N:ns+12+3*N);
yhais      = param(ns+13+3*N:ns+12+4*N);
lis        = param(ns+13+4*N:ns+12+5*N);
mcis       = param(ns+13+5*N:ns+12+6*N);
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);
ls         = sum(lis)/N;
ws         = -UL/UC;
ps         = mus/MPs;

alph1      = 1-alpha;
delt1      = 1-delta;
ome2       = 1-ome1;
rho1       = 1-rho;
thet1      = 1-theta;
mvec       = mus.^(theta/thet1*[0:N-1]);
phbs       = ps/(ome1^(1/rho1)+ome2^(1/rho1))^(-rho1/rho);
o          = ones(N,1);
o1         = ones(N+1,1);

%______________________________________________________________________________
%
% VARIABLES APPEARING IN THE RESIDUALS
%

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

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

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

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

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

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

c          = (ys*y-sum( (xi.*(o*xis'))' )'/N)/cs;
ca         = (ys*ya-sum( (xai.*(o*xis'))' )'/N)/cs;
           % c  = [c_t,...,c_{t+N-1}]'

phi        = zeros(N);
pfi        = zeros(N);
pfai       = zeros(N);
phai       = zeros(N);
mui        = zeros(N);
muai       = zeros(N);
for i=1:N;
  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)'; 
  mui(i,:)   = [0;cumsum(mu(N+i-2:-1:i))]';
  muai(i,:)  = [0;cumsum(mua(N+i-2:-1:i))]';
end;
           % phi = [ph_{t-1},  ph_{t-2},   ..., ph_{t-N},
           %        ph_t    ,  ph_{t-1},   ..., ph_{t-N+1},
           %           :            :             :
           %        ph_{t+N-2},ph_{t+N-3}, ..., ph_{t-1}]
           %
           % mui = [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+N-2},    mu_{t+N-2}+mu_{t+N-3}+...+mu_t]

phb        = (phi-mui)*mvec'/sum(mvec);
pfb        = (pfi-mui)*mvec'/sum(mvec);
pfab       = (pfai-muai)*mvec'/sum(mvec);
phab       = (phai-muai)*mvec'/sum(mvec);
           % phb = [phb_{t-1}, ..., phb_{t+N-2}]'

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+N-1}]'

yhi        = (y+p/rho1-(rho-theta)*phb/rho1/thet1)*o'-(phi-mui)/thet1;
yfi        = (y+p/rho1-(rho-theta)*pfb/rho1/thet1)*o'-(pfi-mui)/thet1;
yfai       = (ya+pa/rho1-(rho-theta)*pfab/rho1/thet1)*o'-(pfai-muai)/thet1;
yhai       = (ya+pa/rho1-(rho-theta)*phab/rho1/thet1)*o'-(phai-muai)/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}]

li         = (yhi.*(o*yhis') + yhai.*(o*yhais'))./(alph1*o*(yhis'+yhais'))- ...
             alpha*ki(1:N,[N,1:N-1])/alph1;
lai        = (yfai.*(o*yhis')+ yfi.*(o*yhais'))./(alph1*o*(yhis'+yhais'))- ...
             alpha*kai(1:N,[N,1:N-1])/alph1;
           % li =  [l_{1,t}      l_{2,t}      ... l_{N,t}
           %        l_{1,t+1}    l_{2,t+1}    ... l_{N,t+1}
           %           :            :                :
           %        l_{1,t+N-1}  l_{2,t+N-1}  ... l_{N,t+N-1}]

l          = li* lis/sum(lis);
la         = lai*lis/sum(lis);
           % l  = [l_t,...,l_{t+N-1}]' 

m          = mu(N:2*N-1)-p;
ma         = mua(N:2*N-1)-pa;
           % m  = [m_t-p_t, ..., m_{t+N-1}-p_{t+N-1}]'

w          = (UCL/UL-UCC/UC)*cs*c+ ...
             (ULL/UL-UCL/UC)*ls*l+  ...
             (ULM/UL-UCM/UC)*MPs*m;
wa         = (UCL/UL-UCC/UC)*cs*ca+ ...
             (ULL/UL-UCL/UC)*ls*la+ ...
             (ULM/UL-UCM/UC)*MPs*ma;
           % w  = [w_t,w_{t+1},...w_{t+N-1}]'

mc         = w+alpha*diag(li)-alpha*diag(ki(1:N,[N,1:N-1]));
mca        = wa+alpha*diag(lai)-alpha*diag(kai(1:N,[N,1:N-1]));
           % mc = [mc_{1,t}, mc_{2,t+1}, ...mc_{N,t+N-1}]'

q          = -UCC*cs/UC*(c-ca)-UCL*ls/UC*(l-la)-UCM*MPs/UC*(m-ma);
           % q  = [q_t,q_{t+1}, ...q_{t+N-1}]'

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

uc         = UCC*cs*c +UCL*ls*l +UCM*MPs*m;
uca        = UCC*cs*ca+UCL*ls*la+UCM*MPs*ma;
           % uc = [U_{c,t},U_{c,t+1}, ... U_{c,t+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);

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

%
% 2. MONOPOLISTS' EULER EQUATIONS FOR CAPITAL
%

tem1           = beta*ws*alpha/alph1*(1-dphis).*lis([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)     =-(UCC*cs /UC*(c(2)-c(1))+ ...
                  UCL*ls /UC*(l(2)-l(1))+ ...
                  UCM*MPs/UC*(m(2)-m(1)) ...
                 +tem1 .* (w(2)+li(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) =-(UCC*cs /UC*(ca(2)-ca(1))+ ...
                  UCL*ls /UC*(la(2)-la(1))+ ...
                  UCM*MPs/UC*(ma(2)-ma(1)) ...
                 +tem1 .* (wa(2)+lai(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/UM-UCC/UC)*cs *c(1)+ ...
                        (ULM/UM-UCL/UC)*ls *l(1)+ ...
                        (UMM/UM-UCM/UC)*MPs*m(1))+ ...
              beta/mus*(UCC*cs /UC*(c(2)-c(1))+ ...
                         UCL*ls /UC*(l(2)-l(1))+ ...
                         UCM*MPs/UC*(m(2)-m(1))+p(1)-p(2)-mu(N));

R(2*N+6,1)     = UM/UC*((UCM/UM-UCC/UC)*cs *ca(1)+ ...
                        (ULM/UM-UCL/UC)*ls *la(1)+ ...
                        (UMM/UM-UCM/UC)*MPs*ma(1))+ ...
              beta/mus*(UCC*cs /UC*(ca(2)-ca(1))+ ...
                         UCL*ls /UC*(la(2)-la(1))+ ...
                         UCM*MPs/UC*(ma(2)-ma(1))+pa(1)-pa(2)-mua(N));

