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

%          Ellen McGrattan, 3-22-00
%          Revised, ERM, 3-22-00


% Subset of parameters used to construct the steady state:
N        = param(2);
Nw       = param(3);
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);
g        = param(36);
v        = param(46);

% Intermediate parameters:
alph1    = 1-alpha;
delt1    = 1-delta;
ome2     = 1-ome1;
thet1    = 1-theta;
v1       = 1-v;
sup1     = 1/(1-rho);
sup2     = rho/(rho-1);
oi       = ones(N,1);
oj       = ones(Nw,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        = .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,UC,UM,UL] ...
         = util(c,MP,l,param);
  res    = [1-beta-UM/UC;
            mc*Fl+UL/UC/v];

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

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

  z0     = z0 - dres\res;

end;
cm       = z0(1);
l        = z0(2);
c        = cl*l-g;
MP       = c/cm;
[U,UC,UM,UL] ...
         = util(c,MP,l,param);
ph       = 1/(MP*pph);
y        = kl^alpha*l/pph;
ki       = kl*l*oi;
xi       = delta*ki;
yhi      = y*(ome1*pph)^sup1 * oi;
yhai     = y*(ome2*pph)^sup1 * 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*UC(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;
  pPh      = (ome1^sup1+ome2^sup1)^(1/sup2);
  Phph     = sum1^(-thet1/theta);
  pph      = pPh*Phph;
  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*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,UC,UM,UL] ...
             = util(cj,mdj,Lsj,param);
    rwage    = -Ww/v *  (UL'*vec2) / (UC'*vec3) ;

    res      = [ mci-rwage *oi./Fli;
                 theta-pph * (mci'*vec1)/sum2;
                 1-beta/mu-UM(1)/UC(1);
                 UC(2:Nw)-UC(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,UC,UM,UL] ...
               = util(cj,mdj,Lsj,param);
      rwage    = -Ww/v *  (UL'*vec2) / (UC'*vec3) ;

      dres(:,j)= ([ mci-rwage * oi./Fli;
                    theta-mci'*[(beta*mu^(1/thet1)).^[0:N-1]'/N]*pph/sum2;
                    1-beta/mu-UM(1)/UC(1);
                    UC(2:Nw)-UC(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,UC,UM,UL] ...
           = util(cj,mdj,Lsj,param);
  rwage    = -Ww/v *  (UL'*vec2) / (UC'*vec3) ;
  p        = mu/mean(mdj);
  wbar     = rwage*p;
  w        = wbar/Ww;
  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(cj,mdj,Lsj,param);
ss       = [ph;y;w;cj;mdj;Lsj;UC;UM;UL;UCC;UMM;ULL;UCM;UCL;ULM;
            ki;xi;yhi;yhai;Ldi;mci];


