%RBCQDE.M
%  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@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.

clear variables

% Specify Parameter Values
% Christiano and Eichenbaum (1992)
%b1 = .9926;     % beta
%b2 = .3440;     % alpha
%b3 = .00374;    % theta
%b4 = .0040;     % gamma
%b7 = .0210;     % delta
%b9 = .0180;     % sigma (gamma)
%b10 = 190.8;    % g
%b11 = .9600;    % rho (g)
%b12 = .0210;    % sigma (g)

% 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] = Gamma * 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) );
Gamma = zeros(2,2);
Gamma(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,dim2] = size(A);
[dim3,dim4] = size(Gamma);

% 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;
while crit1 >= eps1 | crit2 >= eps2
   Fi = (eye(dim1)-B*C)\B;
   Ci = (eye(dim1)-B*C)\A;
   crit1 = max(max(abs(Fi-F))); crit2 = max(max(abs(Ci-C)));
   C = Ci; F = Fi;
   iter = iter+1;
   if iter > 500, 
      disp(' The brute-force iterative procedure did not converge after ')
      disp(' 500 iterations. See Binder and Pesaran (1995, 1997) for alternative ')
      disp(' algorithms to compute the matrix C. '),
   end
end

% 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);
while i <= N
   fp1 = F^i/(eye(dim1)-B*C)/Chat*D1*Gamma^i;
   fp2 = F^i/(eye(dim1)-B*C)/Chat*D2*Gamma^(i+1);
   aux3a = fp1+aux3a;
   aux3b = fp2+aux3b;
   i = i+1;
end
crit3 = max(max(abs(fp1+fp2)));

while crit3 > eps3
   N = N+1;
   fp1 = F^N/(eye(dim1)-B*C)/Chat*D1*Gamma^N;
   fp2 = F^N/(eye(dim1)-B*C)/Chat*D2*Gamma^(N+1);
   aux3a = fp1+aux3a;
   aux3b = fp2+aux3b;
   crit3 = max(max(abs(fp1+fp2)));
end
H = aux3a+aux3b;

% Display Results
disp('The decision rule is: '),
disp('x[t] = C*x[t-1] + H*w[t], '),
disp('where: '), C, H
disp(' and the vector x consists of capital, hours, ')
disp(' consumption, output, and investment, ')
disp(' and w contains the technology shock ') 
disp(' and government expenditure. ')
disp(' The forecast horizon, N, is equal to: '), N
