@  SGBOWDEN.PRG                                                         @
@  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 "at" 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.     @

new;
output file=sgbowden.out reset;

@ 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] = R * 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);
R = 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 = inv(Chat)*Bhat;
A = inv(Chat)*Ahat;

dim1 = rows(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,dim1);                  @ M[0] @
M[dim1+1:2*dim1,.] = eye(dim1);                  @ M[1] @
i = 2;
do 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)*inv(B);
   i = i+1;
endo;
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] @
do while crit3 > eps3;
   M = M|zeros(dim1,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)*inv(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;
   do while i <= N+1;
      K[i*dim1+1:(i+1)*dim1,.] = inv(B)*M[(N+2-i)*dim1+1:(N+3-i)*dim1,.]*B;
      i = i+1;
   endo;
   i = 1;
   do while i <= N+1;
      Y[i*dim1+1:(i+1)*dim1,.] = K[(N+1)*dim1+1:(N+2)*dim1,.]
                                 *inv(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,.]*inv(Chat)
                                     *(D1*matpow(R,N+1-i)+D2*matpow(R,N+2-i));
      i = i+1;
   endo;
   crit3 = maxc(maxc(abs(Y[dim1+1:2*dim1,.]*inv(Chat)*
           (D1*matpow(R,N)+D2*matpow(R,N+1)))));
   N = N+1;
endo;
N = N-1;

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

@ Display Results @
format 9,4;
"The decision rule is:";
"x[t] = C*x[t-1] + H*w[t],";
"where:";?;"C = ";?;;C;?;"H = ";?;;H;?;
"and the vector x consists of consumption and capital,";
"and w contains the technology shock.";
"The forecast horizon, N, is equal to:";?;"N = ";?;?;;N;