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

clear variables

% Specify Parameter Values
alpha=-0.5;
eta=0.1;
rho=2;
beta=0.96; 
delta=0.36;

%%%%%%%%%%%%%%%%%%%%%% Model as in Dynare %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% model(linear);
%   y_k1 = e_yk1;
%   y_k2 = e_yk2;
%   y_l1 = e_yl1;
%   y_l2 = e_yl2;
%   (rho)*(c1(+1)-c2(+1))=(rho-eta)*(c1-c2);
%   w=1/beta*w(-1)+(delta*y_k1+(1-delta)*y_l1)-c1 + alpha*(y_k1-y_k2-(z1(-1)-z2(-1)));
%   c1+c2=delta*y_k1+(1-delta)*y_l1 + delta*y_k2+(1-delta)*y_l2;
%   z1=(rho-eta)*c1-rho*c1(+1)+y_k1(+1);
%   z2=(rho-eta)*c2-rho*c2(+1)+y_k2(+1);
% end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% 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]

Chat=eye(5,5)
Chat(1,1)=rho-eta
Chat(1,4)=-(rho-eta)
Chat(2,1)=-(rho-eta)
Chat(3,1)=1
Chat(4,1)=1
Chat(5,4)=-(rho-eta)

Ahat=zeros(5,5)
Ahat(3,2:4)=[-alpha,1/beta,alpha]

Bhat=zeros(5,5)
Bhat(1,1)=rho
Bhat(1,4)=-rho
Bhat(2,1)=-rho
Bhat(4,1)=1
Bhat(4,4)=1
Bhat(5,4)=-rho

D1=zeros(5,4)
D1(3,1:4)=[alpha+delta,1-delta,-alpha,0]
D1(4,1:4)=[delta,1-delta,delta,1-delta]

D2=zeros(5,4)
D2(2,1)=1
D2(5,3)=1

Gamma=zeros(4,4)

% 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);

% Carry Out Recursions
Q = eye(dim1);
R = Chat\(D1*Gamma^N+D2*Gamma^(N+1));
j = 1;
while j <= N
   R = B/Q*R+Chat\(D1*Gamma^(N-j)+D2*Gamma^(N+1-j));
   Q = eye(dim1)-B/Q*A;
   j = j+1;
end

crit3 = 1;
eps3 = 10^(-6);
QN = Q;
RN = R;

while crit3 > eps3
   Q = QN;
   R = RN;
   QN = eye(dim1);
   RN = Chat\(D1*Gamma^(N+1)+D2*Gamma^(N+2));
   j = 1;
   while j <= (N+1)
      RN = B/QN*RN+Chat\(D1*Gamma^(N+1-j)+D2*Gamma^(N+2-j));
      QN = eye(dim1)-B/QN*A;
      j = j+1;
   end
   crit3 = max(max(abs(QN\RN-Q\R)));
   N = N+1;
end

C = QN\A;
H = QN\R;

% 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