@  SGLDU.PRG                                                            @
@  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 "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=sgldu.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): 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);
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 @
captheta = eye(dim1);
capgamma = inv(Chat)*D1*matpow(R,N)+inv(Chat)*D2*matpow(R,N+1);
i = 2;
do while i <= (N+1);
   captheta = eye(dim1)-B*inv(captheta)*A;
   capgamma = inv(captheta)*(inv(Chat)*D1*matpow(R,N+1-i)
              +inv(Chat)*D2*matpow(R,N+2-i)+B*capgamma);
   i = i+1;
endo;

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

do while crit3 > eps3;
   capgamma = capgammaN;
   capthetaN = eye(dim1);
   capgammaN = inv(Chat)*D1*matpow(R,N+1)+inv(Chat)*D2*matpow(R,N+2);
   i = 2;
   do while i <= (N+2);
      capthetaN = eye(dim1)-B*inv(capthetaN)*A;
      capgammaN = inv(capthetaN)*(inv(Chat)*D1*matpow(R,N+2-i)
                  +inv(Chat)*D2*matpow(R,N+3-i)+B*capgammaN);
      i = i+1;
   endo;
   crit3 = maxc(maxc(abs(capgammaN-capgamma)));
   N = N+1;
endo;

C = inv(capthetaN)*A;
H = capgammaN;

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