function foc=intang_resid(X,param);
%Intang_resid  Residuals for the 3-sector model with intangible capital.
%

%              Ellen R. McGrattan, 3-18-05
%              Revised, ERM, 8-10-05

%______________________________________________________________________________
%
% PARAMETERS

beta       = param(1);
deltI      = param(2);
deltT      = param(3);
eta        = param(4);
gamma      = param(5);
phi1       = param(6);
phi2       = param(7);
psi        = param(8);
sigma      = param(9);
thet1      = param(10);
thet2      = param(11);
kT0        = param(12);
kI0        = param(13);
T          = param(14);
inputs     = zeros(T,15);
inputs(:)  = param(15:14+15*T);
z1         = inputs(:,1);
z2         = inputs(:,2);
xwedge     = inputs(:,3);
g          = inputs(:,4);
tauc       = inputs(:,5);
gtauh      = inputs(:,6);
gtaux      = inputs(:,7);
tauk       = inputs(:,8);
taup       = inputs(:,9);
taud       = inputs(:,10);
taug       = inputs(:,11);
chi        = inputs(:,12);
yn         = inputs(:,13);
xn         = inputs(:,14);
hn         = inputs(:,15);
bderiv     = param(15+15*T);
uderiv     = param(16+15*T);

bhat       = beta*(1+gamma)^(-sigma);
lshare1    = 1-thet1-phi1;
lshare2    = 1-thet2-phi2;
grate      = (1+eta)*(1+gamma);
taux       = (1+gtaux)./xwedge-1;

%______________________________________________________________________________
%
% VARIABLES APPEARING IN RESIDUALS

x          = zeros(T,4);
x(:)       = X;
kTn        = x(:,1);
kIn        = x(:,2);
h1         = x(:,3);
h2         = x(:,4);
kT         = [kT0;kTn(1:T-1)];
kI         = [kI0;kIn(1:T-1)];
xT         = grate*kTn-(1-deltT)*kT;
xI         = grate*kIn-(1-deltI)*kI;
h          = h1+h2+hn;
qxIyb      = lshare1/lshare2 *h2./h1;
k1T        = thet1*kT./(thet1+thet2*qxIyb);
k2T        = kT-k1T;
yb         = k1T.^thet1.*kI.^phi1.*(z1.*h1).^lshare1;
qxI        = qxIyb.*yb;
q          = qxI./xI;
c          = yb+yn-xT-xn-g;
l          = 1-h;
uc         = c.^(-sigma).*l.^(psi*(1-sigma));
rT         = thet1*yb./k1T;
rI         = (phi1*yb+phi2*qxI)./kI;
t          = 1:T-1;
t1         = 2:T;
xtau       = (1+taux).*(1-taud)./(1+tauc);
dtau       = (1-taud).*taup./(1+tauc);
ptau       = (1-taup).*(1-taud)./(1+tauc);
gtau       = q.*(taug+(1-chi).*(1-gtauh)+chi.*(1-taup).*(1-taud))./(1+tauc);

foc        = [[xtau(t).*uc(t)- ...
                   bhat*uc(t1).*(ptau(t1).*(rT(t1)-tauk(t1))+ ...
                            (1-deltT)*xtau(t1)+deltT*dtau(t1));
                kTn(T)-bderiv*kTn(T-1)];
              [gtau(t).*uc(t)- ...
                   bhat*uc(t1).*((1-deltI)*gtau(t1)+ptau(t1).*rI(t1));
                kIn(T)-uderiv*kIn(T-1)];
              psi*(1+tauc).*c.*h1./yb-lshare1*(1-gtauh).*l; 
              xI-k2T.^thet2.*kI.^phi2.*(z2.*h2).^lshare2];

