function ss = steady(param);
%STEADY    calculates the steady states for the models in Chari,
%          Kehoe, and McGrattan, ``Can Sticky Price Models Generate
%          Volatile and Persistent Real Exchange Rates?''
%          Preferences are given in UTIL.m
%

%          Ellen McGrattan, 9-16-99
%          Revised, ERM, 1-30-00


% Subset of parameters used to construct the steady state:
N        = param(2);
Utype    = param(4);
beta     = param(11);
ome1     = param(12);
rho      = param(13);
theta    = param(14);
alpha    = param(15);
delta    = param(16);
b        = param(17);
mu       = param(19);

% Intermediate parameters:
alph1    = 1-alpha;
delt1    = 1-delta;
ome2     = 1-ome1;
thet1    = 1-theta;
sup1     = 1/(1-rho);
sup2     = rho/(rho-1);
o        = ones(N,1);

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

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

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

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

  z0     = z0 - dres\res;


end;
cm       = z0(1);
l        = z0(2);
c        = cl*l;
MP       = c/cm;
ph       = 1/(MP*pph);
y        = kl^alpha*l/pph;
ki       = kl*l*o;
xi       = delta*ki;
yhi      = y*(ome1*pph)^sup1 * o;
yhai     = y*(ome2*pph)^sup1 * o;
li       = l*o;
mci      = mc*o;

% Now turn on mu:
if (mu > 1);
  sum1     = sum((mu^(theta/thet1)).^[0:N-1]')/N;
  sum2     = sum((beta*mu^(theta/thet1)).^[0:N-1]')/N;
  pPh      = (ome1^sup1+ome2^sup1)^(1/sup2);
  Phph     = sum1^(-thet1/theta);
  pph      = pPh*Phph;

  z0       = [ki;y;cm];
  del      = diag(abs(z0))*.0001;
  dres     = zeros(N+2);

  for i=1:100;

    ki     = z0(1:N);
    y      = z0(N+1);
    cm     = z0(N+2);

    kil    = ki([N,1:N-1]);
    Fi     = y*pph*(Phph)^(theta/thet1) * mu.^([0:N-1]'/thet1);
    li     = (Fi./kil.^alpha).^(1/alph1);
    l      = sum(li)/N;
    Fki    = alpha*Fi./kil;
    Fli    = alph1*Fi./li;
    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    = (o./(1-dphi)-beta*(delt1-phip+dphip.*xi./kil)./(1-dphip))./ ...
             (beta*Fki);
    c      = y-sum(xi)/N;
    MP     = c/cm;
    [U,UC,UM,UL] ...
           = util(c,MP,l,param);
    res    = [ mci+UL/UC *o./Fli;
               theta-mci'*[(beta*mu^(1/thet1)).^[0:N-1]'/N]*pph/sum2;
               1-beta/mu-UM/UC];
    for j=1:N+2;
      zp   = z0 + del(:,j);
      ki   = zp(1:N);
      y    = zp(N+1);
      cm   = zp(N+2);
    
      kil  = ki([N,1:N-1]);
      Fi   = y*pph*(Phph)^(theta/thet1) * mu.^([0:N-1]'/thet1);
      li   = (Fi./kil.^alpha).^(1/alph1);
      l    = sum(li)/N;
      Fki  = alpha*Fi./kil;
      Fli  = alph1*Fi./li;
      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  = (o./(1-dphi)-beta*(delt1-phip+dphip.*xi./kil)./(1-dphip))./ ...
             (beta*Fki);
      c    = y-sum(xi)/N;
      MP   = c/cm;
      [U,UC,UM,UL] ...
           = util(c,MP,l,param);
      dres(:,j) ....
           = ([ mci+UL/UC * o./Fli;
                theta-mci'*[(beta*mu^(1/thet1)).^[0:N-1]'/N]*pph/sum2;
                1-beta/mu-UM/UC]-res)/ del(j,j);
    end;

    z0     = z0 - dres\res;

  end;
  ki       = z0(1:N);
  y        = z0(N+1);
  cm       = z0(N+2);

  kil      = ki([N,1:N-1]);
  Fi       = y*pph*(Phph)^(theta/thet1) * mu.^([0:N-1]'/thet1);
  li       = (Fi./kil.^alpha).^(1/alph1);
  l        = sum(li)/N;
  Fki      = alpha*Fi./kil;
  Fli      = alph1*Fi./li;
  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      = (o./(1-dphi)-beta*(delt1-phip+dphip.*xi./kil)./(1-dphip))./ ...
             (beta*Fki);
  c        = y-sum(xi)/N;
  MP       = c/cm;
  p        = mu/MP;
  ph       = p/pph;
  yhi      = y*(ome1*pPh)^sup1*(Phph)^(1/thet1) * mu.^([0:N-1]'/thet1);
  yhai     = y*(ome2*pPh)^sup1*(Phph)^(1/thet1) * mu.^([0:N-1]'/thet1);

end;


[U,UC,UM,UL,UCC,UMM,ULL,UCM,UCL,ULM] ...
         = util(c,MP,l,param);
ss       = [ph;y;c;MP;UC;UM;UL;UCC;UMM;ULL;UCM;UCL;ULM;ki;xi;yhi;yhai;li;mci];

