@  RBCQDE.PRG                                                            @
@  This program solves the Christiano and Eichenbaum (1992)              @
@  RBC-model using the quadratic determinantal                           @
@  equation method of Binder and Pesaran (1995, 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 papers 'Multivariate Linear Rational Expectations Models and      @
@  Macroeconomic Modelling: A Review and Some New Results' and           @
@  '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=rbcqde.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);
dim3 = rows(R);

@ Compute Matrix C Using Brute-Force Iterative Procedure            @
C = eye(dim1);          @ Initial Condition                         @
F = eye(dim1);          @ Initial Condition                         @
eps1 = 10^(-6);         @ Convergence Criterion for F               @
eps2 = 10^(-6);         @ Convergence Criterion for C               @
crit1 = 1; crit2 = 1;   @ Initial Conditions                        @
iter = 0;
do while crit1 >= eps1 or crit2 >= eps2;
   Fi = inv(eye(dim1)-B*C)*B;
   Ci = inv(eye(dim1)-B*C)*A;
   crit1 = maxc(maxc(abs(Fi-F)));
   crit2 = maxc(maxc(abs(Ci-C)));
   C = Ci;
   F = Fi;
   iter = iter+1;
   if iter > 500;
      "The brute-force iterative procedure did not converge after";
      "500 iterations. See Binder and Pesaran (1995, 1997) for alternative";
      "algorithms to compute the matrix C.";
   endif;
endo;

@ Use Recursive Method of Binder and Pesaran (1995) to Compute the    @
@ Forward Part of the Solution - Determine N                          @
eps3 = 10^(-6);    @ Convergence Criterion @
i = 0;
aux3a = zeros(dim1,dim3);
aux3b = zeros(dim1,dim3);
do while i <= N;
   fp1 = matpow(F,i)*inv(eye(dim1)-B*C)*inv(Chat)
         *D1*matpow(R,i);
   fp2 = matpow(F,i)*inv(eye(dim1)-B*C)*inv(Chat)
         *D2*matpow(R,i+1);
   aux3a = fp1+aux3a;
   aux3b = fp2+aux3b;
   i = i+1;
endo;
crit3 = maxc(maxc(abs(fp1+fp2)));

do while crit3 > eps3;
   N = N+1;
   fp1 = matpow(F,N)*inv(eye(dim1)-B*C)*inv(Chat)
         *D1*matpow(R,N);
   fp2 = matpow(F,N)*inv(eye(dim1)-B*C)*inv(Chat)
         *D2*matpow(R,(N+1));
   aux3a = fp1+aux3a;
   aux3b = fp2+aux3b;
   crit3 = maxc(maxc(abs(fp1+fp2)));
endo;
H = aux3a+aux3b;
@ 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;