function [rule,sim] = qr2lq(param, method, nsim, nobs);
%QR2LQ This file accompanies the article ``Solving the Stochastic Growth
%      Model with a Finite Element Method''
%      rule = qr2lq(param,method)  uses three methods to compute 
%      an approximate solution to:
%     
%        max  sum beta^t (exp(w[t])*k[t-1]^alpha-i[t])^(1-tau)/(1-tau)
%      {i[t]}  t               
%
%        subject to  k[t]     = (1-delta)*k[t-1]+i[t]
%                    ome[t+1] = rho*ome[t] + eps[t+1], 
%                    E[eps[t]]= 0, E[eps[t]eps[t]']=sigma
%
%      which has solution, i[t] = h(k[t-1],1,omega[t]).  The
%      program uses a second-order approximation of the objective 
%      function and then computes the optimal linear rule, e.g.  
%
%        i[t]= a0*k[t]+a1+a2*omega[t],  
%
%      The inputs are the parameter vector, param=[beta,alpha,tau,rho,
%      sigma,delta], and  method, e.g.
%          method  = 1         if direct iteration is used
%                  = 2         if a doubling algorithm is used
%                  = 3         if Vaughan's method is used
%      If [rule,sim] = qr2lq(param,method,nsim,nobs) is used instead,
%      then nsim simulations of length nobs are run; sim(:,1) are means
%      and sim(:,2) are standard deviations for k[t],ome[t],i[t].

%      Ellen McGrattan, 6-14-93

% Parameters:

alpha = param(1);
beta  = param(2);
delta = param(3);
rho   = param(4);
sigma = param(5);
tau   = param(6);


% State vector              x[t]   = (k[t-1],1,w[t])'
% Decision variables        u[t]   = i[t]
% Law of motion for x       x[t+1] = A * x[t] + B * u[t]

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

% Steady state of x         k, 1, 0 

k = (alpha*beta/(1-beta*(1-delta)))^(1/(1-alpha));


% Objective function        U  = (exp(x3)*x2^alpha-u)^(1-tau)/(1-tau) 
% Approximation               ~= x'*Q*x+u'*R*u+2*x'*W*u,   x=[x1,x2,x3]'

upar= [alpha,tau];
z   = [k,1,0,delta*k]';
del = diag([k*1e-4,1e-4,1e-6,delta*k*1e-4+1e-6]);
for i=1:4;
  du(i,1) = .5*(qr2u(z+del(:,i),upar)-qr2u(z-del(:,i),upar))/del(i,i);
  for j=1:i;
    du2(i,j) = .25*( qr2u(z+del(:,i)+del(:,j),upar)- ...
                     qr2u(z+del(:,i)-del(:,j),upar)- ...
                     qr2u(z-del(:,i)+del(:,j),upar)+ ...
                     qr2u(z-del(:,i)-del(:,j),upar))/del(i,i)/del(j,j);
  end;
end;
du2 = du2+tril(du2,-1)';
e   = zeros(4,1); 
e(2)= 1;
M   = e*(qr2u(z,upar)-du'*z+.5*z'*du2*z)*e'+ ...
                    .5*(e*du'+du*e'-e*z'*du2-du2*z*e'+du2);
Q = M(1:3,1:3);
W = M(1:3,4);
R = M(4,4);


% 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  = 3;

% 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;

disp(' ')
disp(' ')
disp(sprintf('Check on solution: delta*kbar = %g, rule*xbar = %g', ...
              delta*k,rule*[k;1;0]))
