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

%                 Ellen R. McGrattan, 2-1-00
%                 Revised, ERM, 2-12-00 

%______________________________________________________________________________
%
% PARAMETERS


N          = param(2);
beta       = param(11);
ome1       = param(12);
rho        = param(13);
theta      = param(14);
alpha      = param(15);
delta      = param(16);
b          = param(17);

alph1      = 1-alpha;
delt1      = 1-delta;
ome2       = 1-ome1;
rho1       = 1-rho;
thet1      = 1-theta;
o          = ones(N,1);
o1         = ones(N+1,1);

%______________________________________________________________________________
%
% VARIABLES APPEARING IN THE RESIDUALS
%

Z          = exp(Z);
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)'];
pf         = [X(2*N-2:-1:N);z(2,N:-1:1)'];
pfa        = [X(3*N-3:-1:2*N-1);z(3,N:-1:1)'];
pha        = [X(4*N-4:-1:3*N-2);z(4,N:-1:1)'];
           % 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)'];
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+N-1},  k_{2,t+N-1}, ...  k_{N,t+N-1}]

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

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

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

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

kk         = ki(2:N+1,:) ./ki(1:N,[N,1:N-1]);
kka        = kai(2:N+1,:)./kai(1:N,[N,1:N-1]);

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

c          = y-sum(xi')'/N-g;
ca         = ya-sum(xai')'/N-ga;
           % 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,:)   = [1;cumprod(mu(N+i-2:-1:i))]';
  muai(i,:)  = [1;cumprod(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 = [1, mu_{t-1}, ..., mu_{t-1}*mu_{t-2}*...*mu_{t-N+1},
           %        1, mu_t,          mu_t    *mu_{t-1}*...*mu_{t-N+2} 
           %           :               :             
           %        1, mu_{t+N-2},    mu_{t+N-2}*mu_{t+N-3}*...*mu_t]

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

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

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

li         = ( (yhi+yhai)./ki(1:N,[N,1:N-1]).^alpha ).^(1/alph1) ./(a*ones(1,N));
lai        = ( (yfai+yfi)./kai(1:N,[N,1:N-1]).^alpha ).^(1/alph1)./(aa*ones(1,N));
           % 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          = sum(li')'/N;
la         = sum(lai')'/N;
           % 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}]'

[U,UC,UM,UL]     = util(c,m,l,param);
[Ua,UCa,UMa,ULa] = util(ca,ma,la,param);

w          = -UL./UC;
wa         = -ULa./UCa;
           % w  = [w_t,w_{t+1},...w_{t+N-1}]'

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

q          = UCa./UC;
           % 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 = [lam_{h,t},lam_{h,t+1}, ...lam_{h,t+N-1}]'
%______________________________________________________________________________
%
% LINEARIZED RESIDUALS:
%
% 1. PRICING EQUATIONS
%

bm1        = beta.^[0:N-1]' .* ([1;cumprod(mu(N:2*N-2))]).^(1/thet1);
bm2        = beta.^[0:N-1]' .* ([1;cumprod(mu(N:2*N-2))]).^(theta/thet1);
bma1       = beta.^[0:N-1]' .* ([1;cumprod(mua(N:2*N-2))]).^(1/thet1);
bma2       = beta.^[0:N-1]' .* ([1;cumprod(mua(N:2*N-2))]).^(theta/thet1);
 
R(1,1)     = ph(N)*theta *sum(UC ./p     .*lh .*bm2) -sum(UC .*mc .*lh .*bm1);
R(2,1)     = pf(N)*theta *sum(UCa./(p.*q).*lf .*bm2) -sum(UCa.*mca.*lf .*bm1);
R(3,1)     = pfa(N)*theta*sum(UCa./pa    .*lfa.*bma2)-sum(UCa.*mca.*lfa.*bma1);
R(4,1)     = pha(N)*theta*sum(UC .*q./pa .*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)*o./(1-dphi) - beta*UC(2)*(w(2)*alpha/alph1* ...
                li(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)*(wa(2)*alpha/alph1* ...
                lai(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)/UC(1)  +beta*UC(2)/UC(1)*p(1)/p(2)/mu(N)-1;
R(2*N+6,1) = UMa(1)/UCa(1)+beta*UCa(2)/UCa(1)*pa(1)/pa(2)/mua(N)-1;
