function ss = galiss(param);
%GALISS     calculates the steady states for the model in GALI.m
%

%           Ellen McGrattan, 3-29-04
%           Revised, ERM, 3-30-04


% Subset of parameters used to construct the steady state:
N        = param(2);
Nw       = param(3);
beta     = param(11);
alpha    = param(12);
delta    = param(13);
b        = param(14);
theta    = param(15);
v        = param(16);
mu       = param(17);
g        = param(18);
inhabit  = 0;

% Intermediate parameters:
alph1    = 1-alpha;
delt1    = 1-delta;
thet1    = 1-theta;
v1       = 1-v;
oi       = ones(N,1);
oj       = ones(Nw,1);

% Start with a good guess (mu=1 case):
mc       = theta;
Fk       = (1-beta*delt1)/(beta*mc);
kl       = (Fk/alpha)^(-1/alph1);
Fl       = alph1*kl^alpha;
cl       = kl^alpha-delta*kl;
cm       = 1;
l        = .25;
z0       = [cm;l];
del      = diag(z0)*.0001;

for i=1:10;
  cm     = z0(1);
  l      = z0(2);
  c      = cl*l-g;
  MP     = c/cm;
  [U,U1,U2,UM,UL] ...
         = util(c,c,MP,l,1,param);
  zet    = U1+beta*U2*inhabit;
  res    = [1-beta-UM/zet;
            mc*Fl+UL/zet/v];

  zp     = z0 + del(:,1);
  cm     = zp(1);
  l      = zp(2);
  c      = cl*l-g;
  MP     = c/cm;
  [U,U1,U2,UM,UL] ...
         = util(c,c,MP,l,1,param);
  zet    = U1+beta*U2*inhabit;
  dres(:,1) ...
         = ([1-beta-UM/zet;
             mc*Fl+UL/zet/v] - res)/del(1,1);

  zp     = z0 + del(:,2);
  cm     = zp(1);
  l      = zp(2);
  c      = cl*l-g;
  MP     = c/cm;
  [U,U1,U2,UM,UL] ...
         = util(c,c,MP,l,1,param);
  zet    = U1+beta*U2*inhabit;
  dres(:,2) ...
         = ([1-beta-UM/zet;
             mc*Fl+UL/zet/v] - res)/del(2,2);

  z0     = z0 - dres\res;

end;
cm       = z0(1);
l        = z0(2);
c        = cl*l-g;
MP       = c/cm;
[U,U1,U2,UM,UL] ...
         = util(c,c,MP,l,1,param);
zet      = U1+beta*U2*inhabit;
p        = 1/MP;
y        = kl^alpha*l;
ki       = kl*l*oi;
xi       = delta*ki;
yi       = y*oi;
Ldi      = l*oi;
mci      = mc*oi;
cj       = c*oj; 
mdj      = MP*oj; 
Lsj      = l*oj;
p        = mu/mean(mdj);
w        = -p*UL(1)/(v*zet(1));

% Now turn on mu:
if (mu > 1);
  sum1     = mean((mu^(theta/thet1)).^[0:N-1]');
  sum2     = mean((beta*mu^(theta/thet1)).^[0:N-1]');
  sum3     = mean((mu^(v/v1)).^[0:Nw-1]');
  vec1     = (beta*mu^(1/thet1)).^[0:N-1]'/N;
  vec2     = (beta*mu^(1/v1)).^[0:Nw-1]'/Nw;
  vec3     = (beta*mu^(v/v1)).^[0:Nw-1]'/Nw;
  Pp       = sum1^(-thet1/theta);
  Ww       = sum3^(-v1/v);
  z0       = [ki;y;cj(2:Nw);mdj(1:Nw)];
  del      = diag(abs(z0))*.0001;
  dres     = zeros(N+2*Nw);

  for i=1:100;

    ki       = z0(1:N);
    y        = z0(N+1);
    cj(2:Nw) = z0(N+2:N+Nw);
    mdj      = z0(N+Nw+1:N+2*Nw);

    kil      = ki([N,1:N-1]);
    Fi       = y*(Pp)^(theta/thet1) * mu.^([0:N-1]'/thet1);
    Ldi      = (Fi./kil.^alpha).^(1/alph1);
    Ld       = mean(Ldi);
    Fki      = alpha*Fi./kil;
    Fli      = alph1*Fi./Ldi;
    xi       = (1+b*delta-sqrt(1+2*b-2*b*ki./kil))/b.*kil;
    xil      = xi([N,1:N-1]);
    kil2     = ki([N-1,N,1:N-2]);
    phip     = b/2*(xi./kil-delta).^2;
    dphi     = b*(xil./kil2-delta);
    dphip    = b*(xi./kil-delta);
    mci      = (oi./(1-dphi)-beta*(delt1-phip+dphip.*xi./kil)./(1-dphip))./ ...
               (beta*Fki);
    cj(1)    = Nw*(y-mean(xi)-g)-sum(cj(2:Nw));
    Lsj      = (Ww*mu.^[0:Nw-1]').^(1/v1) * Ld;
    
    [U,U1,U2,UM,UL] ...
             = util(cj,cj,mdj,Lsj,ones(Nw,1),param);
    zet      = U1+beta*U2*inhabit;
    rwage    = -Ww/v *  (UL'*vec2) / (zet'*vec3) ;

    res      = [ mci-rwage *oi./Fli;
                 theta-(mci'*vec1)/sum2;
                 1-beta/mu-UM(1)/zet(1);
                 U1(2:Nw)-U1(1);
                 UM(2:Nw)-UM(1)];

    for j=1:N+2*Nw;
      zp       = z0 + del(:,j);
      ki       = zp(1:N);
      y        = zp(N+1);
      cj(2:Nw) = zp(N+2:N+Nw);
      mdj      = zp(N+Nw+1:N+2*Nw);
    
      kil      = ki([N,1:N-1]);
      Fi       = y*pph*(Phph)^(theta/thet1) * mu.^([0:N-1]'/thet1);
      Ldi      = (Fi./kil.^alpha).^(1/alph1);
      Ld       = mean(Ldi);
      Fki      = alpha*Fi./kil;
      Fli      = alph1*Fi./Ldi;
      xi       = (1+b*delta-sqrt(1+2*b-2*b*ki./kil))/b.*kil;
      xil      = xi([N,1:N-1]);
      kil2     = ki([N-1,N,1:N-2]);
      phip     = b/2*(xi./kil-delta).^2;
      dphi     = b*(xil./kil2-delta);
      dphip    = b*(xi./kil-delta);
      mci      = (oi./(1-dphi)-beta*(delt1-phip+dphip.*xi./kil)./(1-dphip))./ ...
                 (beta*Fki);
      cj(1)    = Nw*(y-mean(xi)-g)-sum(cj(2:Nw));
      Lsj      = (Ww*mu.^[0:Nw-1]').^(1/v1) * Ld;
    
      [U,U1,U2,UM,UL] ...
               = util(cj,cj,mdj,Lsj,ones(Nw,1),param);
      zet      = U1+beta*U2*inhabit;
      rwage    = -Ww/v *  (UL'*vec2) / (zet'*vec3) ;

      dres(:,j)= ([ mci-rwage * oi./Fli;
                    theta-mci'*[(beta*mu^(1/thet1)).^[0:N-1]'/N]*pph/sum2;
                    1-beta/mu-UM(1)/zet(1);
                    U1(2:Nw)-zet(1);
                    UM(2:Nw)-UM(1)]-res)/ del(j,j);
    end;

    z0     = z0 - dres\res;

  end;
  ki       = z0(1:N);
  y        = z0(N+1);
  cj(2:Nw) = z0(N+2:N+Nw);
  mdj      = z0(N+Nw+1:N+2*Nw);

  kil      = ki([N,1:N-1]);
  Fi       = y*pph*(Phph)^(theta/thet1) * mu.^([0:N-1]'/thet1);
  Ldi      = (Fi./kil.^alpha).^(1/alph1);
  Ld       = mean(Ldi);
  Fki      = alpha*Fi./kil;
  Fli      = alph1*Fi./Ldi;
  xi       = (1+b*delta-sqrt(1+2*b-2*b*ki./kil))/b.*kil;
  xil      = xi([N,1:N-1]);
  kil2     = ki([N-1,N,1:N-2]);
  phip     = b/2*(xi./kil-delta).^2;
  dphi     = b*(xil./kil2-delta);
  dphip    = b*(xi./kil-delta);
  mci      = (oi./(1-dphi)-beta*(delt1-phip+dphip.*xi./kil)./(1-dphip))./ ...
             (beta*Fki);
  cj(1)    = Nw*(y-mean(xi)-g)-sum(cj(2:Nw));
  Lsj      = (Ww*mu.^[0:Nw-1]').^(1/v1) * Ld;
    
  [U,U1,U2,UM,UL] ...
           = util(cj,cj,mdj,Lsj,ones(Nw,1),param);
  zet      = U1+beta*U2*inhabit;
  rwage    = -Ww/v *  (UL'*vec2) / (zet'*vec3) ;
  p        = mu/mean(mdj);
  wbar     = rwage*p;
  w        = wbar/Ww;
  yi       = y*(Pp)^(1/thet1) * mu.^([0:N-1]'/thet1);

end;

ss       = [p;w;ki;cj;mdj;y];


