%SGLDU.M
%  This program solves the stochastic growth model
%  in Binder and Pesaran (1998) using the LDU factorization 
%  procedure as described in Proposition 3.2 of Binder and 
%  Pesaran (1999).
%
%  If you find any error in this program, please send e-mail to: 
%  binder@glue.umd.edu
%  Copyright: Michael Binder and M. Hashem Pesaran
%  Current Version: 03/02/2000
%
%  The paper 'Solution of Multivariate Linear Rational Expectations
%  Models and Large Sparse Linear Systems' describing the theory 
%  underlying this program can be downloaded from:
%  http://www.inform.umd.edu/econ/mbinder
%  Please see the user notes at this URL before using this program.

clear variables 

% Specify Parameter Values
% Binder and Pesaran (1998)
b1 = 1.03^(-.25);     % beta
b2 = .8;              % eta
b3 = .64;             % alpha
b4 = .02;             % delta
b5 = .15;             % rho0
b6 = .96;             % rho1
b7 = .02;             % sigma

N = 300;  % Initial Forecasting Horizon (See Binder and Pesaran, 1998)

% Non-Stochastic Steady States
AA = exp(b5/(1-b6)+.5*b7^2/(1-b6^2));
k = ((b1*(1-b3)*AA)/(1-b1*(1-b4)))^(1/b3);
c = AA*k^(1-b3)-b4*k;
sc = c/(c+b4*k);
sk = k/(c+b4*k);

% Set Up the System
% Chat * x[t] = Ahat * x[t-1] + Bhat * E(x[t+1]|I[t])
%               + D1 * w[t] + D2 * E(w[t+1]|I[t])
% w[t] = Gamma * w[t-1] + v[t]

% Variables are in Logs, and in Deviations From Steady State
% First Row of x(t): Consumption
% Second Row of x(t): Capital
% First Row of w(t): Technology Shock

% Coefficient Matrices
Chat = zeros(2,2);
Ahat = zeros(2,2);
Bhat = zeros(2,2);
Chat(1,1) = sc;
Chat(1,2) = sk;
Chat(2,1) = b2;
Chat(2,2) = -b3*b1*(1-b3)*AA*k^(-b3);
Ahat(1,2) = (1-b4)*sk+(1-b3);
Bhat(2,1) = b2;
D1 = zeros(2,1);
D2 = zeros(2,1);
D1(1,1) = 1;
D2(2,1) = -b1*AA*(1-b3)*k^(-b3);
Gamma = b6;

% Transform System to Canonical Form:
% x[t] = A * x[t-1] + B * E(x[t+1]|I[t]) 
%        + inv(Chat) * D1 * w[t] + inv(Chat) * D2 * E(w[t+1]|I[t])
B = Chat\Bhat;
A = Chat\Ahat;

[dim1,dim2] = size(A);

% Carry Out Recursions
captheta = eye(dim1);
capgamma = Chat\D1*Gamma^N+Chat\D2*Gamma^(N+1);
i = 2;
while i <= (N+1)
   captheta = eye(dim1)-B/captheta*A;
   capgamma = captheta\(Chat\D1*Gamma^(N+1-i) ...
              +Chat\D2*Gamma^(N+2-i)+B*capgamma);
   i = i+1;
end

crit3 = 1;
eps3 = 10^(-6);
capgammaN = capgamma;

while crit3 > eps3
   capgamma = capgammaN;
   capthetaN = eye(dim1);
   capgammaN = Chat\D1*Gamma^(N+1)+Chat\D2*Gamma^(N+2);
   i = 2;
   while i <= (N+2)
      capthetaN = eye(dim1)-B/capthetaN*A;
      capgammaN = capthetaN\(Chat\D1*Gamma^(N+2-i) ...
                  +Chat\D2*Gamma^(N+3-i)+B*capgammaN);
      i = i+1;
   end
   crit3 = max(max(abs(capgammaN-capgamma)));
   N = N+1;
end

C = capthetaN\A;
H = capgammaN;

% Display Results
disp('The decision rule is: '),
disp('x[t] = C*x[t-1] + H*w[t], '),
disp('where: '), C, H
disp(' the vector x consists of consumption and capital, ')
disp(' and w contains the technology shock. ')
disp (' The forecast horizon, N, is equal to: '), N
