function [res,V] = transfocc(X,param)
%TRANFOC Computes first-order conditions for balanced growth steady state 
%        for the following model
%
%        max  sum    beta^t (log c[i,t] + psi log (1-l[i,t])) N[i,t]
%           t=0:inf
%
%        subject to
%
%          C[i,t] + Xk[i,t] + Xm[i,t] + NX[i,t]  < = Y[i,t]
%          Y[i,t]    = A[i,t] [N[i,t] hatM[i,t]]^phi
%                            *[K[i,t]^alpha L[i,t]^(1-alpha)]^{1-phi}
%          K[i,t+1]  = (1-deltk) K[i,t] + Xk[i,t]
%          M[i,t+1]  = (1-deltm) M[i,t] + Xm[i,t]
%          hatM[i,t] = M[i,t] + ome[i,t] M[-i,t]
%          NX[i,t]   = sum_{j\=i} (r[i,j,t] M[j,t] - r[j,i,t] M[i,t])
%          M[i,t]   >= 0 
%          N[i,t]    = N[i,0]*(1+gamN)^t
%          A[i,t]    = (1+gamA)^t,  gamA=func(gamY,gamN)
%
%        and in equilibrium,
%          r[j,i,t]  = phi*ome[j,t] Y[j,t]/hatM[j,t]
%
%        Parameters:  I       number of countries 
%                     T       number of periods
%                     alpha   capital share for non-technology input
%                     beta    discount factor
%                     deltk   depreciation rate of non-technology capital
%                     deltm   dpereciation rate of technology capital
%                     gamN    growth rate of N[i,t], all i
%                     gamY    growth rate of Y[i,t], all i
%                     phi     share of technology capital in production
%                     psi     utility parameter
%                     N       Ix1 vector of initial populations
%                     A       T*I x 1 vector of TFPs
%                     theta   T*I x 1 vector of openness parameters 
%

%   Ellen McGrattan, 3-14-07

I          = param(1);
alpha      = param(2);
beta       = param(3);
deltk      = param(4);
deltm      = param(5);
gamN       = param(6);
gamY       = param(7);
phi        = param(8);
psi        = param(9);
N          = param(11:I+10);
ki0        = param(I+11:2*I+10);
mi0        = param(2*I+11:3*I+10);
T          = param(3*I+11);
z          = zeros(T,I);
A          = z;
theta      = z;
A(:)       = param(3*I+12:3*I+11+I*T);
theta(:)   = param(3*I+12+I*T:3*I+11+2*I*T);
nik        = param(3*I+12+2*I*T);   
nim        = param(3*I+13+2*I*T);   
ikt        = param(3*I+14+2*I*T:3*I+13+2*I*T+nik);          % i=find(xk[t]>0)
imt        = param(3*I+14+2*I*T+nik:3*I+13+2*I*T+nik+nim);  % i=find(xm[t]>0)
gamy       = (1+gamY)/(1+gamN)-1;
beth       = beta/(1+gamy);
ome        = theta.^(1/phi);


%
% Unknowns: c[i,t], ydivhm[i,t], l[i,0]
%

c          = z;
l          = z;
xk         = z;
xm         = z;
c(:)       = X(1:I*T);
l(:)       = X(I*T+1:2*I*T);
xk(ikt)    = X(2*I*T+1:2*I*T+nik);
xm(imt)    = X(2*I*T+nik+1:2*I*T+nik+nim);

%
% Intermediate variables
%
k          = z;
k(1,:)     = ki0';
m          = z;
m(1,:)     = mi0';
for i=2:T;
  k(i,:)   = ((1-deltk)*k(i-1,:)+xk(i-1,:))/(1+gamY);
  m(i,:)   = ((1-deltm)*m(i-1,:)+xm(i-1,:))/(1+gamY);
end;
oi         = ones(I,1);
for t=1:T;
  Ome1       = oi*ome(t,:)-diag(ome(t,:))+eye(I);
  Ome2       = ome(t,:)'*oi'-diag(ome(t,:))+eye(I);
  mhat(t,:)  = m(t,:)*Ome2';
  y(t,:)     = A(t,:).*mhat(t,:).^phi.*k(t,:).^(alpha*(1-phi)).* ...
               l(t,:).^((1-alpha)*(1-phi));
  rhs(t,:)   = (N'.*y(t,:)./mhat(t,:))*Ome1';
  r          = ome(t,:).*N'.*y(t,:)./mhat(t,:);
  for i=1:I;
    noti     = [1:i-1,i:I]';
    nx(t,i)  = ome(t,i)*N(i)*y(t,i)/mhat(t,i)* ...
               sum(m(t,noti))-m(t,i)*sum(r(noti));
  end;
end;
nx         = phi*nx;

%
% Residuals
%

N          = ones(T,1)*N';
o          = ones(T-1,I);
res        = [y-psi*c.*l./((1-alpha)*(1-phi)*(1-l)), ...
              N.*c+N.*xk+xm+nx-N.*y];
kres       = [o-beth*c(1:T-1,:)./c(2:T,:).* ...
                (1-deltk+alpha*(1-phi)*y(2:T,:)./k(2:T,:));
              xk(T-1,:)-xk(T,:)];
mres       = [o-beth*c(1:T-1,:)./c(2:T,:).*(1-deltm+phi*rhs(2:T,:));
              m(T-1,:)-m(T,:)];
res        = [res(:);kres(ikt);mres(imt)];

V          = [y(:),c(:),xk(:),xm(:),l(:),k(:),mhat(:),m(:),nx(:)];

