% THE OPTIMAL LINEAR REGULATOR
% The purpose of this program is to clarify and explain the 
% linear quadratic method as described in LJUNGQVIST and SARGENT (CHAPTER 4)
% This program computes the value function and the optimal
% decision rules of the linear-quadratic approximation
% THE ORIGINAL PROGRAM IS WRITTEN BY Jorge Duran e-mail: xurxo@eco.uc3m.es
% THIS PROGRAM IS AMENDED TO INCLUDE THE "USING THE CONSTANT METHOD" AS
% IN SARGENT BOOK
% ADDITION WRITTEN BY O. Mikhail (2002)
% e-mail: omikhail@cleopatra.bus.ucf.edu // omikhail@hotmail.com

%  This message is for the loop below.
diary lqproblem.txt
disp(' The program is iterating on Bellman`s equation. Please wait...')

%  Step 0:  Input the parameter values
  a = 0.33               %  alpha
  b = 0.96               %  beta
  p = 0.95               %  rho
  d = 0.10               %  delta
  e = exp(1)             %  Number e (intrascendent but necessary)

%  Step 1:  Compute the steady state
% The state vector is [K X Z]

  K = ((a*b)/(1-b*(1-d)))^(1/(1-a))
  X = d*K
  Z = 0.0
  
%  Step 2:  Construct the quadratic expansion of the utility function
%  Step 2.1: Evaluate all the first and second order derivatives of the utiliy
%            function at the steady state:
  
% This is the utility function

  R = log(e^Z*K^a - X)

% Now derive (Jacobian)

  Jz = (e^Z*K^a)/(e^Z*K^a - X)
  Jk = (e^Z*a*K^(a-1))/(e^Z*K^a - X)
  Jx = (-1)/(e^Z*K^a - X)

% Get the second derivatives (Hessian)

  Hzz = (((e^Z*K^a)*(e^Z*K^a - X)) - (e^Z*K^a)^2 )/((e^Z*K^a - X)^2)
  Hkk = (((e^Z*a*(a-1)*K^(a-2))*(e^Z*K^a - X)) - ...
                             (e^Z*a*K^(a-1))^2 )/((e^Z*K^a - X)^2)
  Hxx = (-1)/((e^Z*K^a - X)^2)
  Hzk = (((e^Z*a*K^(a-1))*(e^Z*K^a - X)) - (e^Z*K^a*e^Z*a*K^(a-1)) )/ ...
                                                      ((e^Z*K^a - X)^2)
  Hzx = (e^Z*K^a)/((e^Z*K^a - X)^2)
  Hkx = (e^Z*a*K^(a-1))/((e^Z*K^a - X)^2)

%  Step 2.2: Define the Jacobian vector and the Hessian matrix evaluated
%            at the steady state
  DJ  = [Jz  Jk  Jx ]'

  D2H = [Hzz Hzk Hzx
         Hzk Hkk Hkx
         Hzx Hkx Hxx]
         
%  Step 2.3: Define matrix Q
% W here is the Z vector in Sargent

  W = [Z K X]'

% Make the M matrix in Sargent

  Q11 = R - W'*DJ + 0.5*W'*D2H*W
  Q12 = 0.5*(DJ-D2H*W)
  Q22 = 0.5*D2H

  Q=[ Q11  Q12'
      Q12  Q22]

%  Step 3:  Compute the optimal value function.
%  Step 3.1:  Partition matrix Q to separate the state and control variables

  Qff = Q(1:3,1:3)
  Qfd = Q(4,1:3)
  Qdd = Q(4,4)

%  Step 3.2:  Input matrix B

  B = [ 1 0  0  0
        0 p  0  0
        0 0 1-d 1]

%  Step 3.3:  Input matrix P0

  P = [-0.1      0      0
         0    -0.1      0
         0       0   -0.1  ]

%  Step 3.4:  Initialize auxiliary matrix A

  A=ones(3)

%  Step 3.5:  Iterate on Bellman's equation until convergence

  while (norm(A-P)/norm(A))>0.0000000001
    A = P;
    M=B'*P*B; % this is the Bellman equation
      Mff = M(1:3,1:3);
      Mfd = M(4,1:3);
      Mdd = M(4,4);
% The RICATTI equation
    P=Qff + (b*Mff) - (Qfd'+(b*Mfd)')*inv(Qdd+(b*Mdd))*(Qfd+(b*Mfd));
  end
% IF you want to look at the matrices
Mff
Qfd'
b
Mfd
Qdd
Mdd
Qfd
%  Step 4:  Output the results
  disp(' ')
  disp(' The linear policy for investment is d*=J''F, where vector J is:')
  J=(inv(Qdd+(b*Mdd))*(Qfd+(b*Mfd)));
  J
  disp(' ')
  disp(' The optimal value function is V*=F''PF, where matrix P is:')
  P
%  END OF SUPPLIED PROGRAM


%  OR you can do it using the constant as in LJUNGQVIST and SARGENT
%  let the games begin

%  NOW YOU CAN DO IT AS IN LJUNGQVIST and SARGENT
%  THIS IS WRITTEN BY O. MIKHAIL
%  ARE U TALKIN' TO ME, ARE U TALKIN' TO ME, ARE U TALKIN' TO ME

% I do not like using abbreviations for the parameters
  beta = b   
  delta = d
  rho = p

% you will notice that all the matrices here have double letters
% like DJJ D2HH AA BB MM vv RR HH QQ PP JJ and so on
% this does not reflect any bad physical vision
% just that we do not get confused with the matrices in the above program

  DJJ  = [0 Jz  Jk  Jx ]'

  D2HH = [0 0   0   0 
         0 Hzz Hzk Hzx
         0 Hzk Hkk Hkx
         0 Hzx Hkx Hxx]

  WW = [1 Z K X]' % vector of State and Control
  vv = [1 0 0 0]'  % vector e as defined in Sargent
  MM1 = vv*(R-(DJJ'*WW)+(0.5*WW'*D2HH*WW))*vv'
  MM2 = 0.5*(DJJ*vv'-vv*WW'*D2HH-D2HH*WW*vv'+vv*DJJ')
  MM3 = 0.5*D2HH
  MM = MM1+MM2+MM3 % the M matrix in Sargent  
  
  % compare the MM matrix with the Q matrix above
  MM
  Q
  % Now partition the MM matrix into the quadratic format
  RR = MM(1:3,1:3)
  HH = MM(4,1:3)
  QQ = MM(4,4)
  
  AA = [ 1 0   0 
         0 rho 0
         0 0   1-delta]
  BB = [0 0 1]'
  
  [PP,JJ] = olrp(beta,AA,BB,QQ,RR,HH)
  % compare between 
  % the "no constant method P " and 
  % the "constant method PP "
  disp(' The constant method PP')
  PP
  disp(' The no constant method P')
  P
  % compare between 
  % the "no constant method optimal policy rule J " and 
  % the "constant method optimal policy rule JJ "
  disp(' The constant method optimal policy rule JJ ')
  JJ
  disp(' The no constant method optimal policy rule J')
  J
  
  disp(' May the force be with you , BYE BYE NOW ')
  disp(' O. Mikhail ')
  diary off