function x1 = ss_fullmodel(x0,param);
%FULLMODEL_SS    Computes steady state for the full model.
%

n           = param(1);
alpha       = param(2);
beta        = param(3);
cbar        = param(4);
deli        = param(5:4+n);
delt        = param(5+n:4+2*n);
g           = param(5+2*n:4+3*n);
gamma       = zeros(n);
gamma(:)    = param(5+3*n:4+n*n+3*n);
gn          = param(5+n*n+3*n);
gz          = param(6+n*n+3*n);
nu          = zeros(n);
nu(:)       = param(7+n*n+3*n:6+2*n*n+3*n);
omega       = param(7+2*n*n+3*n:6+2*n*n+4*n);
p1          = param(7+2*n*n+4*n);
phi         = param(8+2*n*n+4*n:7+2*n*n+5*n);
psi         = param(8+2*n*n+5*n);
rho         = param(9+2*n*n+5*n);
tauc        = param(10+2*n*n+5*n);
taud        = param(11+2*n*n+5*n);
tauh        = param(12+2*n*n+5*n);
tauk        = param(13+2*n*n+5*n);
taup        = param(14+2*n*n+5*n);
taux        = param(15+2*n*n+5*n);
theta       = param(16+2*n*n+5*n:15+2*n*n+6*n);
z1          = param(16+2*n*n+6*n:15+2*n*n+7*n);
z2          = param(16+2*n*n+7*n:15+2*n*n+8*n);
zeta        = zeros(n);
zeta(:)     = param(16+2*n*n+8*n:15+3*n*n+8*n);
o           = ones(n,1);
betg        = beta*(1+gz)^(-alpha);
gam         = sum(gamma)';
lshare      = 1-theta-phi-gam;
thet1       = 1-theta;
del         = max(abs(x0*1e-4),1e-8);

for i=1:10;
  aggcnrm   = x0(1);
  pvec      = [p1;x0(2:n)];
  h2h1      = x0(n+1:2*n);
  h         = x0(2*n+1:3*n);

  pmvec     = gam.*prod(((pvec*o')./(gamma+1e-15)).^(gamma./(o*gam')))';
  ptvec     = prod(((pvec*o')./(zeta+1e-15)).^zeta)';
  aggp      = sum(omega.^rho.*pvec.^(1-rho))^(1/(1-rho));
  cnrm      = (aggp*omega./pvec).^rho * aggcnrm;
  c         = cnrm*cbar;
  aggc      = (sum(omega.*c.^((rho-1)/rho)))^(rho/(rho-1));
  aggl      = 1-sum(h);
  m1y       = gam.*pvec./pmvec;
  h1        = h./(1+h2h1);
  h2        = h-h1;
  qxipy     = h2h1;
  xiy       = qxipy.^(1-phi).*(z2./z1).^lshare;
  qp        = qxipy./xiy;
  qvec      = qp.*pvec;
  pivec     = prod(((qvec*o')./(nu+1e-15)).^nu)';
  kiy       = betg*phi.*(qvec.*xiy+pvec)./((1-betg*(1-deli)).*pivec);
  ret       = ((1+taux)*ptvec./pvec-betg*(1-(1-taup)*(delt+tauk) + ...
               (1-delt).*((1+taux)*ptvec./pvec-1)))/(betg*(1-taup));
  k1ty      = theta./ret;
  y         = (k1ty.^theta.*kiy.^phi.*m1y.^gam).^(o./lshare).*z1.*h1;
  k1t       = k1ty.*y;
  k2t       = h2h1.*k1t;
  kt        = k1t+k2t;
  m1        = m1y.*y;
  xi        = xiy.*y;
  ki        = kiy.*y;
  qxip      = qvec.*xi./pvec;
  prodxt    = ((1+gz)*(1+gn)-1+delt).*kt;
  prodxi    = ((1+gz)*(1+gn)-1+deli).*ki;
  m1mat     = gamma.*((o./pvec)*pvec').*(o*y');
  m2mat     = gamma.*((o./pvec)*pvec').*(o*qxip');
  xtmat     = zeta.*((o./pvec)*ptvec').*(o*prodxt');
  ximat     = nu.*((o./qvec)*pivec').*(o*prodxi');
  aggl      = 1-sum(h);

  res       = [y-sum(xtmat')'-sum(m1mat')'-sum(m2mat')'-g-c;
               xi-sum(ximat')';
               (1-tauh)*aggl*omega.*lshare.*y.*c.^(-1/rho)./h1- ...
                 (1+tauc)*psi*aggc^((rho-1)/rho)];

  for s=1:3*n;
    x1        = x0;
    x1(s)     = x0(s)+del(s);
    aggcnrm   = x1(1);
    pvec      = [p1;x1(2:n)];
    h2h1      = x1(n+1:2*n);
    h         = x1(2*n+1:3*n);

    pmvec     = gam.*prod(((pvec*o')./(gamma+1e-15)).^(gamma./(o*gam')))';
    ptvec     = prod(((pvec*o')./(zeta+1e-15)).^zeta)';
    aggp      = sum(omega.^rho.*pvec.^(1-rho))^(1/(1-rho));
    cnrm      = (aggp*omega./pvec).^rho * aggcnrm;
    c         = cnrm*cbar;
    aggc      = (sum(omega.*c.^((rho-1)/rho)))^(rho/(rho-1));
    aggl      = 1-sum(h);
    m1y       = gam.*pvec./pmvec;
    h1        = h./(1+h2h1);
    h2        = h-h1;
    qxipy     = h2h1;
    xiy       = qxipy.^(1-phi).*(z2./z1).^lshare;
    qp        = qxipy./xiy;
    qvec      = qp.*pvec;
    pivec     = prod(((qvec*o')./(nu+1e-15)).^nu)';
    kiy       = betg*phi.*(qvec.*xiy+pvec)./((1-betg*(1-deli)).*pivec);
    ret       = ((1+taux)*ptvec./pvec-betg*(1-(1-taup)*(delt+tauk) + ...
                 (1-delt).*((1+taux)*ptvec./pvec-1)))/(betg*(1-taup));
    k1ty      = theta./ret;
    y         = (k1ty.^theta.*kiy.^phi.*m1y.^gam).^(o./lshare).*z1.*h1;
    k1t       = k1ty.*y;
    k2t       = h2h1.*k1t;
    kt        = k1t+k2t;
    m1        = m1y.*y;
    xi        = xiy.*y;
    ki        = kiy.*y;
    qxip      = qvec.*xi./pvec;
    prodxt    = ((1+gz)*(1+gn)-1+delt).*kt;
    prodxi    = ((1+gz)*(1+gn)-1+deli).*ki;
    m1mat     = gamma.*((o./pvec)*pvec').*(o*y');
    m2mat     = gamma.*((o./pvec)*pvec').*(o*qxip');
    xtmat     = zeta.*((o./pvec)*ptvec').*(o*prodxt');
    ximat     = nu.*((o./qvec)*pivec').*(o*prodxi');
    aggl      = 1-sum(h);

    dres(:,s) = ([y-sum(xtmat')'-sum(m1mat')'-sum(m2mat')'-g-c;
                  xi-sum(ximat')';
                  (1-tauh)*aggl*omega.*lshare.*y.*c.^(-1/rho)./h1- ...
                    (1+tauc)*psi*aggc^((rho-1)/rho)] - res)/del(s);
  end;
  x0          = x0-real(dres\res);
end;
x1          = x0;


