function [res,V] = bgfocc(X,param)
%BGFOCC Computes first-order conditions for balanced growth steady state 
%       for the following model (with M=0 imposed):
%
%       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  (zeta is penalty parameter)
%         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 
%                   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       Ix1 vector of TFPs
%                   sigma   Ix1 vector of openness parameters (=ome^phi)
%

%   Ellen McGrattan, 3-13-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);
zeta    = param(10);
N       = param(11:I+10);
A       = param(I+11:2*I+10);
sigma   = param(2*I+11:3*I+10);
ni      = param(3*I+11);
imt     = param(3*I+12:3*I+11+ni);
gamy    = (1+gamY)/(1+gamN)-1;
irate   = (1+gamy)/beta-1;
ome     = sigma.^(1/phi);

%
% Unknowns: c, l[0]
%

l       = X(1:I);
m       = zeros(I,1);
m(imt)  = X(I+1:I+ni); 

%
% Intermediate variables
%
o       = ones(I,1);
Ome1    = o*ome'-diag(ome)+eye(I);
Ome2    = ome*o'-diag(ome)+eye(I);
ydivk   = (irate+deltk)/(alpha*(1-phi))*o;
mhat    = Ome2*m;
y       = (A.*mhat.^phi.*ydivk.^(-alpha*(1-phi)).*l.^ ...
          ((1-alpha)*(1-phi))).^(1/(1-alpha*(1-phi)));
ydivhm  = y./mhat;
k       = y./ydivk;
c       = (1-alpha)*(1-phi)*(1-l).*y./l/psi;
xm      = (gamY+deltm)*m;
xk      = (gamY+deltk)*k;
nx      = zeros(I,1);
r       = ome.*N.*ydivhm;
for i=1:I;
  noti  = [1:i-1,i:I]';
  nx(i) = ome(i)*N(i)*y(i)/mhat(i)*sum(m(noti))-m(i)*sum(r(noti));
end;
nx      = phi*nx;

%
% Residuals
%
mres    = 1-beta/(1+gamy)*(1-deltm+phi*Ome1*(N.*ydivhm));
res     = [N.*(y-c-xk)-xm-nx;
           mres(imt)];
V       = [y,c,xk,xm,l,k,mhat,m,nx];
