function R = res_niscmtr_ms(Z,param)
%RES_NISCMTR_MS   Residuals for the international sticky price model with
%                 Taylor rule for interest rates and N>1.

%                 Ellen R. McGrattan, 2-1-00
%                 Revised, ERM, 7-23-01

%______________________________________________________________________________
%
% 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);
mu         = param(19);
irule      = param(23:32);
ys         = param(ns+1);
rs         = mu/beta-1;

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
%

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

ph         = [X(N-1:-1:1);z(1,N:-1:1)'];
pf         = [X(2*N+1:-1:N+3);z(2,N:-1:1)'];
pfa        = [X(3*N+3:-1:2*N+5);z(3,N:-1:1)'];
pha        = [X(4*N+5:-1:3*N+7);z(4,N:-1:1)'];
phl        = [X(N+2),X(N+1),X(N)]';
pfl        = [X(2*N+4),X(2*N+3),X(2*N+2)]';
pfal       = [X(3*N+6),X(3*N+5),X(3*N+4)]';
phal       = [X(4*N+8),X(4*N+7),X(4*N+6)]';
           % ph  = [ph_{t-N}, ...,ph_t, ... ph_{t+N-2}]'   (2*N-1 x 1)
           % phl = [ph_{t-N-3},ph_{t-N-2},ph_{t-N-1}]'
           % similar for pf,pf*,ph*

ki         = [X(4*N+9:5*N+8)';z(5:N+4,N:-1:1)'];
kai        = [X(5*N+9:6*N+8)';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)';

yl1        = X(6*N+9);
yl2        = X(6*N+10);
yal1       = X(6*N+11);
yal2       = X(6*N+12);
           % y   = [y_t, ... y_{t+N-1}]'
           % yli = y_{t-i}

r          = z(2*N+7,N:-1:1)';
ra         = z(2*N+8,N:-1:1)';
rl1        = X(6*N+13);
rl2        = X(6*N+14);
rl3        = X(6*N+15);
ral1       = X(6*N+16);
ral2       = X(6*N+17);
ral3       = X(6*N+18);
           % r   = [r_t, ... r_{t+N-1}]'
           % rli = r_{t-i}

eps        = Z(nz*N+nX+1);
epsa       = Z(nz*N+nX+2);
           % eps = eps_t 

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

a          = exp(Z(nz*N+nX+3*N+2:-1:nz*N+nX+2*N+3));
aa         = exp(Z(nz*N+nX+4*N+2:-1:nz*N+nX+3*N+3));
           % 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);
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)';
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}]
           %

phil1      = [phi(1,2:N),phl(3)];
phil2      = [phil1(2:N),phl(2)];
phil3      = [phil2(2:N),phl(1)];
pfil1      = [pfi(1,2:N),pfl(3)];
pfil2      = [pfil1(2:N),pfl(2)];
pfil3      = [pfil2(2:N),pfl(1)];
pfail1     = [pfai(1,2:N),pfal(3)];
pfail2     = [pfail1(2:N),pfal(2)];
pfail3     = [pfail2(2:N),pfal(1)];
phail1     = [phai(1,2:N),phal(3)];
phail2     = [phail1(2:N),phal(2)];
phail3     = [phail2(2:N),phal(1)];
           % phili = [ph_{t-i-1}, ..., ph_{t-N-i}]

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

phbl1      = (sum( (phil1' ./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
phbl2      = (sum( (phil2' ./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
phbl3      = (sum( (phil3' ./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
pfbl1      = (sum( (pfil1' ./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
pfbl2      = (sum( (pfil2' ./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
pfbl3      = (sum( (pfil3' ./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
pfabl1     = (sum( (pfail1'./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
pfabl2     = (sum( (pfail2'./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
pfabl3     = (sum( (pfail3'./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
phabl1     = (sum( (phail1'./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
phabl2     = (sum( (phail2'./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
phabl3     = (sum( (phail3'./mu1).^(-theta/thet1) )/N)^(-thet1/theta);
           % phbli = phb_{t-i-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);

pl1        = (ome1^(1/rho1)*phbl1  ^(-rho/rho1) ...
             +ome2^(1/rho1)*pfbl1  ^(-rho/rho1))^(-rho1/rho);
pl2        = (ome1^(1/rho1)*phbl2  ^(-rho/rho1) ...
             +ome2^(1/rho1)*pfbl2  ^(-rho/rho1))^(-rho1/rho);
pl3        = (ome1^(1/rho1)*phbl3  ^(-rho/rho1) ...
             +ome2^(1/rho1)*pfbl3  ^(-rho/rho1))^(-rho1/rho);
pal1       = (ome1^(1/rho1)*pfabl1 ^(-rho/rho1) ...
             +ome2^(1/rho1)*phabl1 ^(-rho/rho1))^(-rho1/rho);
pal2       = (ome1^(1/rho1)*pfabl2 ^(-rho/rho1) ...
             +ome2^(1/rho1)*phabl2 ^(-rho/rho1))^(-rho1/rho);
pal3       = (ome1^(1/rho1)*pfabl3 ^(-rho/rho1) ...
             +ome2^(1/rho1)*phabl3 ^(-rho/rho1))^(-rho1/rho);
           % p   = [p_t,p_{t+1}, ..., p_{t+N-1}]'
           % pli =  p_{t-i}

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') .* (mui ./pfai).^(1/thet1);
yhai       = (ya.*(ome2*pa).^(1/rho1).*phab.^((theta-rho)/rho1/thet1)* ...
              o') .* (mui ./phai).^(1/thet1);

yhail1     = (yal1*(ome2*pal1)^(1/rho1)*phabl1^((theta-rho)/rho1/thet1)* ...
              o') .* (mu1'./phail1).^(1/thet1);
yhail2     = (yal2*(ome2*pal2)^(1/rho1)*phabl2^((theta-rho)/rho1/thet1)* ...
              o') .* (mu1'./phail2).^(1/thet1);
yfil1      = (yl1 *(ome2*pl1) ^(1/rho1)*pfbl1 ^((theta-rho)/rho1/thet1)* ...
              o') .* (mu1'./pfil1).^(1/thet1);
yfil2      = (yl2 *(ome2*pl2) ^(1/rho1)*pfbl2 ^((theta-rho)/rho1/thet1)* ...
              o') .* (mu1'./pfil2).^(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}]
           % yfil1 = [yf_{1,t-1}  yf_{2,t-1} ... yf_{N,t-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*o./p;
ma         = mu*o./pa;
z0         = [m;ma];
del        = diag(z0)*.0001;
for i=1:10;
  m            = z0(1:N);
  ma           = z0(N+1:2*N);
  [U,UC,UM]    = util(c,m,l,param);
  [Ua,UCa,UMa] = util(ca,ma,la,param);
  res          = [UM./UC-r./(o+r);
                  UMa./UCa-ra./(o+ra)];

  for j=1:2*N;
    zp         = z0 + del(:,j);
    m          = zp(1:N);
    ma         = zp(N+1:2*N);
    [U,UC,UM]  = util(c,m,l,param);
    [Ua,UCa,UMa] ...
               = util(ca,ma,la,param);
    dres(:,j) ....
               = ([UM./UC-r./(o+r);
                   UMa./UCa-ra./(o+ra)]-res)/ del(j,j);
  end;
  z0           = z0 - dres\res;

end;
m          = z0(1:N);
ma         = z0(N+1:2*N);
           % 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]' .* mu1.^(1/thet1);
bm2        = beta.^[0:N-1]' .* mu1.^(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.*bm2) -sum(UCa.*mca.*lfa.*bm1);
R(4,1)     = pha(N)*theta*sum(UC .*q./pa .*lha.*bm2) -sum(UC .*mc .*lha.*bm1);

phs        = param(ns);        % Normalize residuals 
MPs        = param(ns+3);
UCs        = param(ns+4);
ps         = mu/MPs;
lhs        = (ome1*ps)^(1/rho1)*phs^((theta-rho)/(rho1*thet1))*ys;
R(1)       = R(1)/(phs*theta*UCs*lhs/ps);
R(2)       = R(2)/(phs*theta*UCs*lhs/ps);
R(3)       = R(3)/(phs*theta*UCs*lhs/ps);
R(4)       = R(4)/(phs*theta*UCs*lhs/ps);

%
% 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-1;
R(2*N+6,1) = UMa(1)/UCa(1)+beta*UCa(2)/UCa(1)*pa(1)/pa(2)/mu-1;

%
% 4. INTEREST RATE EQUATIONS
%
gdp          = y(1) + sum(yhai(1,:))/N-sum(yfi(1,:))/N;
gdpl1        = yl1  + sum(yhail1)/N-sum(yfil1)/N;
gdpl2        = yl2  + sum(yhail2)/N-sum(yfil2)/N;
dy           = log(gdp)   -log(ys);
dyl1         = log(gdpl1) -log(ys);
dyl2         = log(gdpl2) -log(ys);
gdpa         = ya(1)+ sum(yfi(1,:))/N-sum(yhai(1,:))/N;
gdpal1       = yal1 + sum(yfil1)/N-sum(yhail1)/N;
gdpal2       = yal2 + sum(yfil2)/N-sum(yhail2)/N;
dya          = log(gdpa)  -log(ys);
dyal1        = log(gdpal1)-log(ys);
dyal2        = log(gdpal2)-log(ys);
pif1         = log(p(2))-log(p(1));
pi           = log(p(1))-log(pl1);
pil1         = log(pl1) -log(pl2);
pil2         = log(pl2) -log(pl3);
piaf1        = log(pa(2))-log(pa(1));
pia          = log(pa(1))-log(pal1);
pial1        = log(pal1) -log(pal2);
pial2        = log(pal2) -log(pal3);

R(2*N+7,1)   = r(1) -rs-irule' *[rl1-rs;rl2-rs;rl3-rs;
                                pif1;pi;pil1;pil2;dy;dyl1;dyl2]-eps;
R(2*N+8,1)   = ra(1)-rs-irule' *[ral1-rs;ral2-rs;ral3-rs;
                                piaf1;pia;pial1;pial2;dya;dyal1;dyal2]-epsa;

