@ RBCFRM.PRG                                                            @
@  This program solves the Christiano and Eichenbaum (1992)             @
@  RBC-model using the fully recursive method of                        @
@  Proposition 4.1 of Binder and Pesaran (1997).                        @
@                                                                       @
@  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 'Multivariate Linear Rational Expectations Models:         @
@  Characterization of the Nature of the Solutions and their Fully      @
@  Recursive Computation' 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=rbcfrm.out reset;

@ Specify Parameter Values @
@ Binder and Pesaran (1995) (Not Rounded) @
b1 = 1.03^(-.25);     @ beta @
b2 = .34474214;       @ alpha @
b3 = .0037398493;     @ theta @
b4 = .0037976621;     @ gamma @
b7 = .020868824;      @ delta @
b9 = .018287903;      @ sigma (gamma) @
rhog0 = .22492956;
rhog1 = .95716082;
sigg = .020638302;
b10 = exp(rhog0/(1-rhog1)+.5*sigg^2/(1-rhog1^2));   @ g @
b11 = rhog1;          @ rho (g) @
b12 = sigg;           @ sigma (g) @

N = 300; @ Initial Forecasting Horizon (See Binder and Pesaran, 1995) @

@ Non-Stochastic Steady States @
g = b10;
sk = b1*b2/(1-b1*(1-b7)*exp(-b4));
H = ((1-b2)/b3+b10/(sk^(b2/(1-b2))*exp(-b2*b4
    /(1-b2))))/(1-(1-(1-b7)*exp(-b4))*sk);
k = (sk*exp(-b2*b4))^(1/(1-b2))*H;
im = k-(1-b7)*exp(-b4)*k;
y = k^b2*H^(1-b2)*exp(-b2*b4);
c = y-im-g;
si = im/y;
sc = c/y;
sg = g/y;

@ 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, 'Detrended', and in Deviations    @
@ From Steady State                                        @
@ First Row of x(t): Capital                               @
@ Second Row of x(t): Hours                                @
@ Third Row of x(t): Consumption                           @
@ Fourth Row of x(t): Output                               @
@ Fifth Row of x(t): Investment                            @
@ First Row of w(t): Technology Shock                      @
@ Second Row of w(t): Government Expenditure               @

@ Coefficient Matrices @
Chat = zeros(5,5);
Ahat = zeros(5,5);
Bhat = zeros(5,5);
Chat[1,3] = -sc;
Chat[1,4] = 1;
Chat[1,5] = -si;
Chat[2,1] = 1;
Chat[2,5] = -(1-exp(-b4)*(1-b7));
Chat[3,2] = -(1-b2);
Chat[3,4] = 1;
Chat[4,2] = b2;
Chat[4,3] = 1;
Chat[5,1] = -b1*b2*(b2-1)*k^(b2-1)*H^(1-b2)*exp(-b2*b4);
Chat[5,3] = -1;
Ahat[2,1] = (1-b7)*exp(-b4);
Ahat[3,1] = b2;
Ahat[4,1] = b2;
Bhat[5,2] = b1*b2*(1-b2)*k^(b2-1)*H^(1-b2)*exp(-b2*b4);
Bhat[5,3] = -1;
D1 = zeros(5,2);
D2 = zeros(5,2);
D1[1,2] = sg;
D1[2,1] = (1-b7)*exp(-b4);
D1[3,1] = b2;
D1[4,1] = b2;
D2[5,1] = b1*( b2^2*k^(b2-1)*H^(1-b2)*exp(-b2*b4)
          +(1-b7)*exp(-b4) );
R = zeros(2,2);
R[2,2] = b11;

@ 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 @
Q = eye(dim1);
RR = inv(Chat)*(D1*matpow(R,N)+D2*matpow(R,N+1));
j = 1;
do while j <= N;
   RR = B*inv(Q)*RR+inv(Chat)*(D1*matpow(R,N-j)+D2*matpow(R,N+1-j));
   Q = eye(dim1)-B*inv(Q)*A;
   j = j+1;
endo;

crit3 = 1;
eps3 = 10^(-6);
QN = Q;
RRN = RR;

do while crit3 > eps3;
   Q = QN;
   RR = RRN;
   QN = eye(dim1);
   RRN = inv(Chat)*(D1*matpow(R,N+1)+D2*matpow(R,N+2));
   j = 1;
   do while j <= (N+1);
      RRN = B*inv(QN)*RRN+inv(Chat)*(D1*matpow(R,N+1-j)
            +D2*matpow(R,N+2-j));
      QN = eye(dim1)-B*inv(QN)*A;
      j = j+1;
   endo;
   crit3 = maxc(maxc(abs(inv(QN)*RRN-inv(Q)*RR)));
   N = N+1;
endo;

C = inv(QN)*A;
H = inv(QN)*RR;

@ 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 capital, hours,";
"consumption, output, and investment,";
"and w contains the technology shock";
"and government expenditure.";
"The forecast horizon, N, is equal to:";?;"N = ";?;?;;N;