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

%              Ellen R. McGrattan, 9-15-99
%              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);
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         = mu/MPs;
rs         = mu/beta-1;

alph1      = 1-alpha;
delt1      = 1-delta;
ome2       = 1-ome1;
rho1       = 1-rho;
thet1      = 1-theta;
mvec       = mu.^(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+8;
nX         = 6*N+18;
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+1:-1:N+3);z(2,N:-1:1)'] - log(phs);
pfa        = [X(3*N+3:-1:2*N+5);z(3,N:-1:1)'] - log(phs);
pha        = [X(4*N+5:-1:3*N+7);z(4,N:-1:1)'] - log(phs);
phl        = [X(N+2),X(N+1),X(N)]' - log(phs);
pfl        = [X(2*N+4),X(2*N+3),X(2*N+2)]' - log(phs);
pfal       = [X(3*N+6),X(3*N+5),X(3*N+4)]' - log(phs);
phal       = [X(4*N+8),X(4*N+7),X(4*N+6)]' - log(phs);
           % 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)'] - log(o1*kis');
kai        = [X(5*N+9:6*N+8)';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);

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

r          = z(2*N+7,N:-1:1)' -rs;
ra         = z(2*N+8,N:-1:1)' -rs;
rl1        = X(6*N+13) -rs;
rl2        = X(6*N+14) -rs;
rl3        = X(6*N+15) -rs;
ral1       = X(6*N+16) -rs;
ral2       = X(6*N+17) -rs;
ral3       = X(6*N+18) -rs;
           % 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 

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);
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}]
           %
phb        = phi*mvec'/sum(mvec);
pfb        = pfi*mvec'/sum(mvec);
pfab       = pfai*mvec'/sum(mvec);
phab       = phai*mvec'/sum(mvec);
           % phb = [phb_{t-1}, ..., phb_{t+N-2}]'

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}]

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

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

yhi        = (y+p/rho1-(rho-theta)*phb/rho1/thet1)*o'-phi/thet1;
yfi        = (y+p/rho1-(rho-theta)*pfb/rho1/thet1)*o'-pfi/thet1;
yfai       = (ya+pa/rho1-(rho-theta)*pfab/rho1/thet1)*o'-pfai/thet1;
yhai       = (ya+pa/rho1-(rho-theta)*phab/rho1/thet1)*o'-phai/thet1;
yhail1     = (yal1+pal1/rho1-(rho-theta)*phabl1/rho1/thet1)*o'-phail1/thet1;
yhail2     = (yal2+pal2/rho1-(rho-theta)*phabl2/rho1/thet1)*o'-phail2/thet1;
yfil1      = (yl1+pl1/rho1-(rho-theta)*pfbl1/rho1/thet1)*o'-pfil1/thet1;
yfil2      = (yl2+pl2/rho1-(rho-theta)*pfbl2/rho1/thet1)*o'-pfil2/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}]
           %
           % yfili = [yf_{1,t-i}  yf_{2,t-i} ... yf_{N,t-i}]

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          = (-(UCM/UM-UCC/UC)*cs*c -(ULM/UM-UCL/UC)*ls*l +beta/mu/rs*r)/ ...
             ((UMM/UM-UCM/UC)*MPs);
ma         = (-(UCM/UM-UCC/UC)*cs*ca-(ULM/UM-UCL/UC)*ls*la+beta/mu/rs*ra)/ ...
             ((UMM/UM-UCM/UC)*MPs);
           % 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*mu^(1/thet1)).^[0:N-1]'.*mcis; w1 = w1/sum(w1);
w2             = (beta*mu^(theta/thet1)).^[0:N-1]';   w2 = w2/sum(w2);

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

%
% 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/mu*(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));

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/mu*(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));

%
% 4. INTEREST RATE EQUATIONS
%


gdp          = y(1) + (yhai(1,:)-yfi(1,:))*yhais/N/ys;
gdpl1        = yl1  + (yhail1-yfil1)*yhais/N/ys;
gdpl2        = yl2  + (yhail2-yfil2)*yhais/N/ys;
gdpa         = ya(1)+ (yfi(1,:)-yhai(1,:))*yhais/N/ys;
gdpal1       = yal1 + (yfil1-yhail1)*yhais/N/ys;
gdpal2       = yal2 + (yfil2-yhail2)*yhais/N/ys;

R(2*N+7,1)   = r(1)-irule' *[rl1;rl2;rl3;p(2)-p(1);p(1)-pl1;
                            pl1-pl2;pl2-pl3;gdp;gdpl1;gdpl2]-eps;
R(2*N+8,1)   = ra(1)-irule'*[ral1;ral2;ral3;pa(2)-pa(1);pa(1)-pal1;
                             pal1-pal2;pal2-pal3;gdpa;gdpal1;gdpal2]-epsa;

