function V1 = balgrowth(N,A,sigma,constant); 
%BALGROWTH.m  Computes the balanced growth steady state for
%
%   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]
%
%   Inputs:  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        = length(N); 
alpha    = constant(1);
beta     = constant(2);
deltk    = constant(3);
deltm    = constant(4);
gamN     = constant(5);
gamY     = constant(6);
phi      = constant(7);
psi      = constant(8);
zeta     = 0;
o        = ones(I,1);
par0     = [I;constant;zeta;N;A;sigma];

%
% Make an educated guess for y/hat m, l
%
ome      = sigma.^(1/phi);
gamy     = (1+gamY)/(1+gamN)-1;
irate    = (1+gamy)/beta-1;
Ome      = o*ome'-diag(ome)+eye(I);
ydivhm   = (Ome\((irate+deltm)/phi*o))./N;
X0       = [1/(1+psi)*o;ydivhm];
X0       = secant('bgfoc',X0,par0);
[r,Vu]   = bgfoc(X0,par0);
m        = Vu(:,8);
imt      = find(m>0);
par0(10) = .01;
for i=1:100;
  par0(10) = par0(10)*1.1;
  X0       = secant('bgfoc',X0,par0);
end;
[r,V]    = bgfoc(X0,par0);
l        = V(:,5);
m        = V(:,8);
X1       = real([l;m(imt)]);
param    = [par0;length(imt);imt];
X1       = secant('bgfocc',X1,param);
[r,V1]   = bgfocc(X1,param);
if any(abs(imag(X1(:)))>0);
  X1     = secant('bgfocc',real(X1),param);
  [r,V1] = bgfocc(X1,param);
end;  
m        = V1(:,8);
if any(m<0); % double check the conditions
  imt    = find(m>0);
  param  = [par0;length(imt);imt];
  X1     = secant('bgfocc',X1,param);
  [r,V1] = bgfocc(X1,param);
end;

