function R = res_niscmtr1_ms(Z,param)
%RES_NISCMTR1_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-25-01

%______________________________________________________________________________
%
% PARAMETERS


ns         = 50;
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;

%______________________________________________________________________________
%
% VARIABLES APPEARING IN THE RESIDUALS
%

nz         = 10;
z          = zeros(10,2);
z(:)       = Z(1:20);
z          = z(:,[2,1]);
z(1:8,:)   = exp(z(1:8,:));
phlag      = exp(Z([23,22,21]));
pflag      = exp(Z([26,25,24]));
pfalag     = exp(Z([29,28,27]));
phalag     = exp(Z([32,31,30]));
           % phlag = [ph_{t-4},ph_{t-3},ph_{t-2}]'

k          = [exp(Z(33));z(5,:)'];
ka         = [exp(Z(34));z(6,:)'];
           % k  = [k_{t-1}, k_t, k_{t+1}]'

ylag       = exp(Z([36,35]));
yalag      = exp(Z([38,37]));
           % ylag  = [y_{t-2},y_{t-1}]'

rlag       = exp(Z([41,40,39]));
ralag      = exp(Z([44,43,42]));
           % rlag  = [r_{t-3},r_{t-2},r_{t-1}]'

eps        = Z(45);
epsa       = Z(46);
           % eps  = eps_t

ph         = z(1,:)';
pf         = z(2,:)';
pfa        = z(3,:)';
pha        = z(4,:)'; 
           % ph  = [ph_{t-1},ph_t]'

y          = z(7,:)'; 
ya         = z(8,:)'; 
           % y   = [y_t,y_{t+1}]'

r          = z(9,:)';
ra         = z(10,:)'; 
           % r   = [r_t,r_{t+1}]'

kk         = k(2:3,:) ./k(1:2);
kka        = ka(2:3,:)./ka(1:2);

x          = (1+b*delta-sqrt(1+2*b-2*b*kk))/b .* k(1:2);
xa         = (1+b*delta-sqrt(1+2*b-2*b*kka))/b.* ka(1:2);
           % x  = [x_t,x_{t+1}]'

c          = y-x;
ca         = ya-xa;
           % c  = [c_t,c_{t+1}]'

p          = (ome1^(1/rho1)*ph .^(-rho/rho1) ...
             +ome2^(1/rho1)*pf .^(-rho/rho1)).^(-rho1/rho);
pa         = (ome1^(1/rho1)*pfa.^(-rho/rho1) ...
             +ome2^(1/rho1)*pha.^(-rho/rho1)).^(-rho1/rho);
           % p  = [p_t,p_{t+1}]'

yh         = y .*(ome1*p ./ph) .^(1/rho1);
yf         = y .*(ome2*p ./pf) .^(1/rho1);
yfa        = ya.*(ome1*pa./pfa).^(1/rho1);
yha        = ya.*(ome2*pa./pha).^(1/rho1);
           % yh = [yh_t,yh_{t+1}]'

l          = ( (yh+yha)./k(1:2) .^alpha ).^(1/alph1) ;
la         = ( (yfa+yf)./ka(1:2).^alpha ).^(1/alph1) ;
           % l  = [l_t,l_{t+1}]'

o          = ones(2,1);
m          = mu*o./p;
ma         = mu*o./pa;
z0         = [m;ma];
del        = diag(z0)*.0001;
for i=1:10;
  m            = z0(1:2);
  ma           = z0(3:4);
  [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:4;
    zp         = z0 + del(:,j);
    m          = zp(1:2);
    ma         = zp(3:4);
    [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:2);
ma         = z0(3:4);
           % m  = [m_t-p_t,m_{t+1}-p_{t+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}]'

mc         = w ./(alph1*(k(1:2) ./l) .^alpha );
mca        = wa./(alph1*(ka(1:2)./la).^alpha );
           % mc = [mc_t,mc_{t+1}]'

q          = UCa./UC;
           % q  = [q_t,q_{t+1}]'

lh         = (ome1*p).^(1/rho1) .*ph .^((theta-rho)/(rho1*thet1)).*y;
lf         = (ome1*p).^(1/rho1) .*pf .^((theta-rho)/(rho1*thet1)).*y;
lfa        = (ome1*pa).^(1/rho1).*pfa.^((theta-rho)/(rho1*thet1)).*ya;
lha        = (ome1*pa).^(1/rho1).*pha.^((theta-rho)/(rho1*thet1)).*ya;
           % lh = [lam_{h,t},lam_{h,t+1}]'

%______________________________________________________________________________
%
% LINEARIZED RESIDUALS:
%
% 1. PRICING EQUATIONS
%


R(1,1)     = ph(1) *theta*UC(1)/p(1)        *lh(1) - UC(1) *mc(1) *lh(1);
R(2,1)     = pf(1) *theta*UCa(1)/(p(1)*q(1))*lf(1) - UCa(1)*mca(1)*lf(1);
R(3,1)     = pfa(1)*theta*UCa(1)/pa(1)      *lfa(1)- UCa(1)*mca(1)*lfa(1);
R(4,1)     = pha(1)*theta*UC(1)*q(1)/pa(1)  *lha(1)- UC(1) *mc(1) *lha(1);

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*(x(2)/k(2)-delta)^2;
dphi       = b*(x(1)/k(1)-delta);
dphin      = b*(x(2)/k(2)-delta);
R(5,1)     = UC(1)/(1-dphi) - beta*UC(2)*(w(2)*alpha/alph1* ...
                l(2)/k(2)+(1-delta-phin+dphin*x(2)/k(2))/(1-dphin));

phin       = b/2*(xa(2)/ka(2)-delta)^2;
dphi       = b*(xa(1)/ka(1)-delta);
dphin      = b*(xa(2)/ka(2)-delta);
R(6,1)     = UCa(1)/(1-dphi) - beta*UCa(2)*(wa(2)*alpha/alph1* ...
                la(2)/ka(2)+(1-delta-phin+dphin*xa(2)/ka(2))/(1-dphin));

%
% 3. MONEY DEMAND EQUATIONS
%

R(7,1)     = UM(1)/UC(1)  +beta*UC(2)/UC(1)*p(1)/p(2)/mu-1;
R(8,1)     = UMa(1)/UCa(1)+beta*UCa(2)/UCa(1)*pa(1)/pa(2)/mu-1;

%
% 4. INTEREST RATE EQUATIONS
%
plag       = (ome1^(1/rho1)*phlag.^(-rho/rho1) ...
             +ome2^(1/rho1)*pflag.^(-rho/rho1)).^(-rho1/rho);
palag      = (ome1^(1/rho1)*pfalag.^(-rho/rho1) ...
             +ome2^(1/rho1)*phalag.^(-rho/rho1)).^(-rho1/rho);
           % plag  = [p_{t-3},p_{t-2},p_{t-1}]'

yflag      = ylag .*(ome2*plag(2:3) ./pflag(2:3)).^(1/rho1);
yhalag     = yalag.*(ome2*palag(2:3)./phalag(2:3)).^(1/rho1);
           % yflag = [yf_{t-2},yf_{t-1}]'

gdp        = y    + yha-yf;
gdpa       = ya   + yf-yha;
gdplag     = ylag + yhalag-yflag;
gdpalag    = yalag+ yflag-yhalag;
           % gdp     = [gdp_t,gdp_{t+1}]'
           % gdplag  = [gdp_{t-2},gdp_{t-1}]'

R(9,1)     = r(1) -rs-irule' *[rlag(3)-rs;rlag(2)-rs;rlag(1)-rs;
                               log(p(2))-log(p(1));
                               log(p(1))-log(plag(3));
                               log(plag(3))-log(plag(2));
                               log(plag(2))-log(plag(1));
                               log(gdp(1))-log(ys);
                               log(gdplag(2))-log(ys);
                               log(gdplag(1))-log(ys)]-eps;

R(10,1)    = ra(1)-rs-irule' *[ralag(3)-rs;ralag(2)-rs;ralag(1)-rs;
                               log(pa(2))-log(pa(1));
                               log(pa(1))-log(palag(3));
                               log(palag(3))-log(palag(2));
                               log(palag(2))-log(palag(1));
                               log(gdpa(1))-log(ys);
                               log(gdpalag(2))-log(ys);
                               log(gdpalag(1))-log(ys)]-epsa;
