function R = res_sticky1(Z,param)
%RES_STICKY1    Residuals for general sticky price model with benchmark
%               functional form for utility (called by sticky1).

%               Ellen R. McGrattan, 1-20-98
%               Revised, ERM, 2-16-99

%______________________________________________________________________________
%
% PARAMETERS

N          = param(1);
psi        = param(2);
sigma      = param(3);
omega      = param(4);
eta        = param(5);
beta       = param(6);
alp1       = param(7);
alp2       = param(8);
delta      = param(9);
b          = param(10);
epsilon    = param(11);
chi        = param(12);

omeg1      = 1-omega;
delt1      = 1-delta;
sigm1      = 1-sigma;
nu         = (eta-1)/eta;
gamma      = 1/(1+psi);
gamm1      = 1-gamma;

%______________________________________________________________________________
%
% VARIABLES APPEARING IN THE RESIDUALS 
%

z          = zeros(N+2,N);
z(:)       = Z(1:(N+2)*N);
plag       = Z((N+2)*N+N-1:-1:(N+2)*N+1);
klag       = Z((N+2)*N+N:(N+2)*N+2*N-1);

pm         = [plag;z(1,N:-1:1)'];
           % pm = [p_{t-N}-m_{t-N}, ...p_{t+N-2}-m_{t+N-2}]'  (2*N-1 x 1)

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

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

Y          = exp(y);
           % Y  = [Y_t, ...Y_{t+N-1}]'

ki         = [klag';z(2:N+1,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}]

KK         = exp(ki(2:N+1,:)-ki(1:N,[N,1:N-1]));
           % KK = [K_{1,t}/K_{N,t-1}, K_{2,t}/K_{1,t-1} ... K_{N,t}/K_{N-1,t-1}
           %      [K_{1,t+1}/K_{N,t}, K_{2,t+1}/K_{1,t} ... K_{N,t+1}/K_{N-1,t}
           %          :            :             :
           %      [K_{1,t+N-1}/K_{N,t+N-2},  ...            K_{N,t+N-1}/
           %                                                     K_{N-1,t+N-2}]
            
Xi         = (1+b*delta-sqrt(1+2*b*delta-2*b*(KK-1+delta)))/b .* ...
             exp(ki(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;
           % C  = [C_t,...C_{t+N-1}]'


nn         = (ones(N,1)*[1:N-1] - triu(N*ones(N,N-1)))/N;
for i=1:N;
  m(i,1)   =-sum(pm(i:i+N-1))/N+[1:N]*mu(i:i+N-1)/N; 
  for j=1:N;
    ppi(i,j) = sum(pm(i:N+i-1))/N-pm(N-j+i)-nn(N-j+1,:)*mu(i:N+i-2);
 end; 
end;
           % m  = [m_t-pb_t,     ... m_{t+N-1}-pb_{t+N-1}]'
           % ppi= [pb_t-p_{t-1}  pb_t-p_{t-2},    ... pb_t-p_{t-N}
           %       pb_{t+1}-p_t  pb_{t+1}-p_{t-1} ... pb_{t+1}-p_{t-N+1}
           %            :                     :                 :
           %       pb_{t+N-1}    pb_{t+N-1}           pb_{t+N-1}
           %         -p_{t+N-2}    -p_{t+N-3} ...       -p_{t-1}  ]

yi         = y*ones(1,N) + epsilon*ppi;
           % yi = [y_{1,t}     y_{2,t}     ... y_{N,t}
           %       y_{1,t+1}   y_{2,t+1}   ... y_{N,t+1}
           %          :            :             :
           %       y_{1,t+N-1} y_{2,t+N-1} ... y_{N,t+N-1}] 

li         = (yi-alp1*ki(1:N,[N,1:N-1]))/alp2;
           % 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(exp(li)')'/N;
           % L  = [L_t,...L_{t+N-1}]'

M          = exp(m);
           % M  = [M_t/Pb_t,     ... M_{t+N-1}/Pb_{t+N-1}]'

U          = ((omega*C.^nu+omeg1*M.^nu).^(gamma/nu).*(1-L).^gamm1).^sigm1/sigm1;
UC         = sigm1*gamma*U*omega.*C.^(nu-1)./(omega*C.^nu+omeg1*M.^nu);
UM         = sigm1*gamma*U*omeg1.*M.^(nu-1)./(omega*C.^nu+omeg1*M.^nu);
UL         = -sigm1*gamm1*U./(1-L);
W          = -UL./UC;
           % W  = [W_t,...W_{t+N-1}]'

mc         = log(W)+(1-alp2)*diag(li)-alp1*diag(ki(1:N,[N,1:N-1]))-log(alp2); 
           % mc = [mc_{1,t},mc_{2,t+1},...mc_{N,t+N-1}]'


XK         = [Xi(1,:)      ./exp(ki(1,[N,1:N-1]));
              Xi(2,[2:N,1])./exp(ki(2,:))];
           % XK = [X_{1,t}/K_{N,t-1}, X_{2,t}/K_{1,t-1},... X_{N,t}/K_{N-1,t-1}
           %       X_{2,t+1}/K_{1,t}, X_{3,t+1}/K_{2,t} ... X_{1,t+1}/K_{N,t}]

phip       = b*(XK(2,:)'-delta).^2/2;
dphi       = -b*delta+b* XK(1,:)';
dphip      = -b*delta+b* XK(2,:)';
           % = adjustment costs and their derivatives

%______________________________________________________________________________
%
% RESIDUALS:
%
% 1. PRICING EQUATION
%

B0         = beta.^[0:N-1]';            % = [1,beta,...,beta^(N-1)]
B1         = cumsum(beta.^[N-1:-1:1]'); % = [beta^(N-1), ...
                                        %    beta^(N-1)+beta^(N-2), ...
                                        %    beta^(N-1)+beta^(N-2)+...+beta]
B2         = cumsum(B1);               
B2         = B2(N-1:-1:1);              % = [beta^(N-1)+
                                        %    beta^(N-1)+beta^(N-2)+...
                                        %    beta^(N-1)+beta^(N-2)+...+beta,
                                        %                  :
                                        %    beta^(N-1)+beta^(N-1)+beta^(N-2),
                                        %    beta^(N-1)]
B1         = B1(N-1:-1:1);              % = [beta+beta^2+...+beta^(N-1),
                                        %         beta^2+...+beta^(N-1), ..., 
                                        %                    beta^(N-1)]
B3         = cumsum(B0(1:N-1));         % = [1,1+beta,...,1+beta+..+beta^(N-2)]
B4         = cumsum(B3);                % = [1,
                                        %    1+(1+beta),...,
                                        %    1+(1+beta)+...+(1+...+beta^(N-2))]
betasum    = sum(B0);                   % = 1+beta+....+beta^(N-1);

varphi     = (1-1/epsilon)/(2-chi/epsilon);
R(1,1)     = (N-1)*betasum*pm(N)-B3'*pm(1:N-1)-B1'*pm(N+1:2*N-1)- ...
             N*varphi*B0'*mc+N*varphi*betasum*log(1-1/epsilon)- ...
             B2'*mu(N:2*N-2)+B4'*mu(1:N-1);          


%
% 2. MONOPOLISTS' EULER EQUATIONS FOR CAPITAL 
%

R(2:N+1,1) = ones(N,1)./(1-dphi)-beta*UC(2)/UC(1)*(W(2)*alp1/alp2* ...
             exp(li(2,[2:N,1])'-ki(2,:)')+ ...
             (delt1-phip+dphip.*XK(2,:)')./(1-dphip));

%
% 3. MONEY DEMAND EQUATION
%

R(N+2,1)   = UM(1)/UC(1)-1+beta*UC(2)/UC(1)*M(2)/M(1)/exp(mu(N+1));

