function z = plant2(z0);
%PLANT2   is a function used to compute the optimal lotteries for
%         the planner in ``The firm and the plant in general equilibrium'' 
%         by Hornstein and Prescott (Section 3).  Here, I assume that
%         all of the mass is on two points.  The problem is:
%
%             max           x*U(c1,h)+(1-x)*U(c0,0)
%           c0,c1,h,x
%                      s.t. x*c1+(1-x)*c0 = h*kbar^alpha*x^(1-alpha)
%
%         z = plant2(z0) takes an initial guess, z0=(c0,h,x) and computes
%         the optimum (z).

%         Ellen McGrattan, 04-06-93



% Inputs
% ------
%
% parameters:

alpha = .3;
gamma = .3;
sigma = .8;
kbar  = 2;



%----------------------------------------------------------------------------

% Problem
% -------

%   max          x*U(c1,h)+(1-x)*U(c0,0)
% c0,c1,h,x
%               s.t. x*c1+(1-x)*c0 = h*k^alpha*x  = h*kbar^alpha*x^(1-alpha)
%
% or
%
%   max          x*U(c1(c0,h,x),h)+(1-x)*U(c0,0)
%  c0,h,x

% First Order Conditions
% ----------------------
%
%  (c0)     x*U_1(c1,h)*c1_1(c0,h,x)+(1-x)*U_1(c0,0) =  0
%  (h)      x*U_1(c1,h)*c1_2(c0,h,x)+x*U_2(c1,h)     =  0
%  (x)      U(c1,h)+x*U_1(c1,h)*c1_3(c0,h,x)-U(c0,0) =  0
%
%  e.g.,  G(y) = 0, y=(c0,h,x) -->  y^{n+1} = y^n - [G'(y^n)]^{-1}G(y^n) 

% Utility
% -------
%
%  Let
%
%  U(c,h)     = ((c^gamma*(1-h)^(1-gamma))^(1-sigma)-1)/(1-sigma)
%  U_1(c,h)   = gamma*c^(gamma*(1-sigma)-1)*(1-h)^((1-gamma)*(1-sigma))
%  U_2(c,h)   = -(1-gamma)*c^(gamma*(1-sigma)-1)*(1-h)^((1-gamma)*(1-sigma))
%  c1(c,h,x)  = h*(kbar/x)^alpha-(1-x)*c/x
%  c1_1(c,h,x)= -(1-x)/x
%  c1_2(c,h,x)= (kbar/x)^alpha
%  c1_3(c,h,x)= -alpha*h*kbar^alpha*x^(-alpha-1)+c/(x*x)

del = .001*ones(3,1);

for i = 1:15;

  c0  = z0(1);
  h   = z0(2);
  x   = z0(3);

  c1 = h*(kbar/x)^alpha-(1-x)*c0/x;
  u0 = (c0^(gamma*(1-sigma)))/(1-sigma);
  u1 = ((c1^gamma*(1-h)^(1-gamma))^(1-sigma))/(1-sigma);
  u11= gamma*(1-sigma)*u1/c1;
  u12= -(1-gamma)*(1-sigma)*u1/(1-h);
  c11= -(1-x)/x;
  c12= (kbar/x)^alpha;
  c13= -alpha*h*kbar^alpha*x^(-alpha-1)+c0/(x*x);
  u01= gamma*(1-sigma)*u0/c0;

  G  = [x*u11*c11+(1-x)*u01;
        x*u11*c12+x*u12;
        u1+x*u11*c13-u0];

  c0 = z0(1)-del(1);

  c1 = h*(kbar/x)^alpha-(1-x)*c0/x;
  u0 = (c0^(gamma*(1-sigma)))/(1-sigma);
  u1 = ((c1^gamma*(1-h)^(1-gamma))^(1-sigma))/(1-sigma);
  u11= gamma*(1-sigma)*u1/c1;
  u12= -(1-gamma)*(1-sigma)*u1/(1-h);
  c11= -(1-x)/x;
  c12= (kbar/x)^alpha;
  c13= -alpha*h*kbar^alpha*x^(-alpha-1)+c0/(x*x);
  u01= gamma*(1-sigma)*u0/c0;

  G1 = [x*u11*c11+(1-x)*u01;
        x*u11*c12+x*u12;
        u1+x*u11*c13-u0];

  c0 = z0(1);
  h  = z0(2)-del(2);

  c1 = h*(kbar/x)^alpha-(1-x)*c0/x;
  u0 = (c0^(gamma*(1-sigma)))/(1-sigma);
  u1 = ((c1^gamma*(1-h)^(1-gamma))^(1-sigma))/(1-sigma);
  u11= gamma*(1-sigma)*u1/c1;
  u12= -(1-gamma)*(1-sigma)*u1/(1-h);
  c11= -(1-x)/x;
  c12= (kbar/x)^alpha;
  c13= -alpha*h*kbar^alpha*x^(-alpha-1)+c0/(x*x);
  u01= gamma*(1-sigma)*u0/c0;

  G2 = [x*u11*c11+(1-x)*u01;
        x*u11*c12+x*u12;
        u1+x*u11*c13-u0];

  h  = z0(2);
  x  = z0(3)-del(3);

  c1 = h*(kbar/x)^alpha-(1-x)*c0/x;
  u0 = (c0^(gamma*(1-sigma)))/(1-sigma);
  u1 = ((c1^gamma*(1-h)^(1-gamma))^(1-sigma))/(1-sigma);
  u11= gamma*(1-sigma)*u1/c1;
  u12= -(1-gamma)*(1-sigma)*u1/(1-h);
  c11= -(1-x)/x;
  c12= (kbar/x)^alpha;
  c13= -alpha*h*kbar^alpha*x^(-alpha-1)+c0/(x*x);
  u01= gamma*(1-sigma)*u0/c0;

  G3 = [x*u11*c11+(1-x)*u01;
        x*u11*c12+x*u12;
        u1+x*u11*c13-u0];

  z0 = z0 - ([(G-G1)/del(1),(G-G2)/del(2),(G-G3)/del(3)])\G;

z0
end;
z=z0;
