alpha = 0.00;
beta  = 0.96;
delta = 0.09;
gamma = 0.25;
lam   = 1.5136;
omega = 1.001;
pa    = 0.7;
pb    = 1.0;
psi   = 0.81;
sig   = 0.021;
theta = 0.25;
zeta  = -.042;
par   = [alpha;beta;delta;gamma;lam;omega;pa;pb;psi;sig;theta;
         0.00003;0.08;0.002];
 

A = [1-delta,0,0,0;
     0,0,0,0;
     0,0,psi,0
     0,0,0,1];
B = [0,1,0,0;
     0,0,1,0;
     0,0,0,0;
     0,0,0,0];

% Steady state: call sscumc.m

ss   = secant('sscumc',[.6;.4;.2;.8],par);


% Objective function Approximation   
%
% U  = (n*u(c1,1-h)+(1-n)*u(c0,1)-p(n)) ~= x'*Q*x+u'*R*u+2*x'*W*u,       
%
% where 
%
%    c0       = (z*f(K,n,h)-cg-i-n*c1-alpha*(n-n_)^2)/(1-n) 
%    u(c,l)   = (c^gamma*l^(1-gamma))^(1-omega)/(1-omega) 
%    f(K,n,h) = K^theta*n^(1-theta)*h
%    x        = [K,n_,ln(z)] 
%    u        = [c1,i,n,h]
%   

z   = [ss(4);ss(2);0;1;ss(3);delta*ss(4);ss(2);ss(1)];
du2 = zeros(8);
del = diag(max(1e-8,abs(z)*1e-4));
for i=1:8;
  du(i,1) = .5*(cumcu(z+del(:,i),par)-cumcu(z-del(:,i),par))/del(i,i);
  for j=1:i;
    du2(i,j) = .25*( cumcu(z+del(:,i)+del(:,j),par)- ...
                     cumcu(z+del(:,i)-del(:,j),par)- ...
                     cumcu(z-del(:,i)+del(:,j),par)+ ...
                     cumcu(z-del(:,i)-del(:,j),par))/del(i,i)/del(j,j);
  end;
end;
du2  = du2+tril(du2,-1)';
e    = zeros(8,1); 
e(4) = 1;
M    = e*(cumcu(z,par)-du'*z+.5*z'*du2*z)*e'+ ...
                    .5*(e*du'+du*e'-e*z'*du2-du2*z*e'+du2);
Q = M(1:4,1:4);
W = M(1:4,5:8);
R = M(5:8,5:8);

method = 1;

% NOTE: to check convergence of the Riccati matrix, follow instructions
%       concerning ``uncommenting'' and ``deleting'' lines.


% Convert to the no-discount, no-cross product case:

rw = R\W';
a  = sqrt(beta)*(A-B*rw);
b  = sqrt(beta)*B;
q  = Q-W*rw;
m  = 4;

% Additional parameters that I set inside the program:
tol   = 1e-8;
maxit = 1000;

if method == 1;                   % Direct Iteration
 
  p0  = -.01*eye(m);
  f0  = (R+b'*p0*b)\(b'*p0*a);
  d1  = 1;
  %d2 = 1;                                       % uncomment
  it  = 1;
  %while ((d1>tol | d2>tol) & it<=maxit);        % uncomment
  while (d1>tol & it<=maxit);                    % delete
    p0  = a'*p0*(a-b*f0) + q;                    % delete
    %p1 = a'*p0*(a-b*f0) + q;                    % uncomment
    f1  = (R+b'*p0*b)\(b'*p0*a);                 % delete
    %f1 = (R+b'*p1*b)\(b'*p1*a);                 % uncomment
    d1  = norm(f1-f0);
    %d2 = norm(p1-p0);                           % uncomment
    it  = it+1;
    f0  = f1;
    %p0 = p1;                                    % uncomment
  end;
  rule = -(f1+rw);
  if it>=maxit; 
    disp(sprintf('WARNING: Iteration limit of %g reached.',maxit)) 
  end;

elseif method == 2                % Doubling Algorithm

  v   = eye(m);
  a0  = a;
  b0  = b*(R\b');
  p0  = q;
  f0  = (b'*p0*b+R)\(b'*p0*a);
  d1  = 1;
  %d2 = 1;
  it  = 1;

  %while ((d1>tol | d2>tol) & it<=maxit);        % uncomment
  while (d1>tol & it<=maxit);                    % delete
    v1  = inv(v+b0*p0);
    a1  = a0*v1*a0;
    b1  = b0+a0*v1*b0*a0';
    p0  = p0+a0'*p0*v1*a0;                       % delete
    %p1 = p0+a0'*p0*v1*a0;                       % uncomment
    f1  = (b'*p0*b+R)\(b'*p0*a);
    %f1 = (b'*p1*b+R)\(b'*p1*a);                 % uncomment
    d1  = norm(f1-f0);
    %d2 = norm(p1-p0);                           % uncomment
    it  = it+1;
    f0  = f1;
    a0  = a1;
    b0  = b1;
    %p0 = p1;                                    % uncomment
  end
  rule = -(f1+rw); 
  if it>=maxit; 
    disp(sprintf('WARNING: Iteration limit of %g reached.',maxit)) 
  end;

else                              % Vaughan's Algorithm

  ai    = inv(a);  
  ri    = inv(R);
  H     = [ai,ai*b*ri*b';
          q*ai,q*ai*b*ri*b'+a'];
  [V,D] = eig(H);

  d    = diag(D);
  i    = find(real(d)==NaN);
  d(i) = inf*ones(length(i),1);
  [d,index] = sort(-abs(d));     % sort on magnitude of eigenvalues
  if (~((-d(m) > 1) & (-d(m+1) < 1)))
      error('Can''t order eigenvalues, (A,B) may be uncontrollable.')
  end
  v11  = V(1:m,index(1:m));
  v21  = V(m+1:2*m,index(1:m));
  p    = real(v21/v11);
  rule = -((R+b'*p*b)\(b'*p*a)+rw);
end;

if nargin > 2;
  
  rand('normal')
  sim = zeros(4,2);
  for i = 1:nsim;
    x   = ltitr(A+B*rule,sigma*C(:,3),rand(nobs,1),[k;1;0]);
    y   = [x,x*rule'];
    sim = sim+[mean(y)',std(y)'];
  end; 
  sim = sim([1,3,4],:)/nsim;

end;

[z(5:8),rule*z(1:4)]
