%SGBOWDEN.M
%  This program solves the stochastic growth model
%  in Binder and Pesaran (1998) using Bowden's procedure
%  as described in Proposition 3.1 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): Capital
% First Row of w(t): Technology Shock

% Coefficient Matrices 
Chat = -b2*sk/sc-b3*b1*(1-b3)*AA*k^(-b3) ...
       -1/sc*((1-b4)*sk+(1-b3))*b2;
Ahat = -b2/sc*((1-b4)*sk+(1-b3));
Bhat = -sk/sc*b2;
D1 = -b2/sc;
D2 = 1/sc*b2-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
eps3 = 10^(-6);
crit3 = 1;
M = zeros((N+3)*dim1,dim1);   % M[0], M[1], ..., M[N+2]
M(1:dim1,:) = zeros(dim1);               % M[0]
M(dim1+1:2*dim1,:) = eye(dim1);          % M[1]
i = 2;
while i <= N+1
   M(i*dim1+1:(i+1)*dim1,:) = (M((i-1)*dim1+1:i*dim1,:) ...
                              -M((i-2)*dim1+1:(i-1)*dim1,:)*A)/B;        
   i = i+1;
end
M((N+2)*dim1+1:(N+3)*dim1,:) = (M((N+1)*dim1+1:(N+2)*dim1,:) ...
                               -M(N*dim1+1:(N+1)*dim1,:)*A);       % M[N+2]
while crit3 > eps3
   M = [M; zeros(dim1)];
   K = zeros((N+3)*dim1,dim1);      % K[0], K[1], ..., K[N+2]
   Y = zeros((N+3)*dim1,dim1);      % Y[0], Y[1], ..., Y[N+2]
   Ystar = zeros((N+3)*dim1,dim1);  % Ystar[0], Ystar[1], ..., 
                                    % Ystar[N+2]
   M((N+1)*dim1+1:(N+2)*dim1,:) = (M(N*dim1+1:(N+1)*dim1,:) ...
                                  -M((N-1)*dim1+1:N*dim1,:)*A)/B;  % M[N+1]
   M((N+2)*dim1+1:(N+3)*dim1,:) = (M((N+1)*dim1+1:(N+2)*dim1,:) ...
                                  -M(N*dim1+1:(N+1)*dim1,:)*A);    % M[N+2]
   K(1:dim1,:) = -M((N+2)*dim1+1:(N+3)*dim1,:);                    % K[0]
   K(dim1+1:2*dim1,:) = M((N+1)*dim1+1:(N+2)*dim1,:);              % K[1]
   i = 2;
   while i <= N+1
      K(i*dim1+1:(i+1)*dim1,:) = B\M((N+2-i)*dim1+1:(N+3-i)*dim1,:)*B;
      i = i+1;
   end
   i = 1;
   while i <= N+1
      Y(i*dim1+1:(i+1)*dim1,:) = K((N+1)*dim1+1:(N+2)*dim1,:) ...
                                 /M((N+2)*dim1+1:(N+3)*dim1,:) ...
                                 *M(i*dim1+1:(i+1)*dim1,:);
      Ystar(i*dim1+1:(i+1)*dim1,:) = Ystar((i-1)*dim1+1:i*dim1,:) ...
                                     +Y(i*dim1+1:(i+1)*dim1,:)/Chat ...
                                     *(D1*Gamma^(N+1-i)+D2*Gamma^(N+2-i)); 
      i = i+1;
   end
   crit3 = max(max(abs(Y(dim1+1:2*dim1,:)/Chat* ...
          (D1*Gamma^N+D2*Gamma^(N+1)))));
   N = N+1;
end
N = N-1;

% Augment the System to Include Consumption
M1 = [1 sk/sc; 0 1];
M2 = [0 1/sc*((1-b4)*sk+(1-b3)); 0 Y((N+1)*dim1+1:(N+2)*dim1,:)*A];
M3 = [1/sc; Ystar((N+1)*dim1+1:(N+2)*dim1,:)]; 
C = M1\M2;
H = M1\M3;

% 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
