%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%            RBC model with fixed labor supply                %
%                                                             %
%                                                             %
%                  For a description see                      %
%              Manuel Wälti (December 2003)                   %
%       "Solving linear rational expectations models"         %
%                                                             %
%            (c) Manuel Waelti, December 2003                 %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                % Clear the memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                     Algorithm Parameters                    %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nfres=1;                  % # of free state variables
npres=2;                  % # of predetermined state variables
ntots=nfres+npres;        % Total # of state variables
ndist=4;                  % # of fundamental disturbances and one-step-ahead deviations
nrep=200;                 % IRF horizon
nshoc=1;                  % # of exogenous shocks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%             Structural Parameters of the Economy            %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta=0.99;                % Discount factor
alpha=0.4;                % Elasticity of output with regard to capital
delta=0.025;              % Depreciation rate
rho=0.95;                 % Autocorrelation of the shock
stda=0.01;                % Standard deviation of the shock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%          Steady-State Values and Parameters a and b         %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
koy=(alpha*beta)/(1-beta*(1-delta)); 
ioy=delta*koy;
coy=1-ioy;
yss=koy^(alpha/(1-alpha)); 
kss=koy*yss;
css=coy*yss;
sss=1;

a1=beta*alpha*(alpha-1)*sss*kss^(alpha-1);
a2=beta*alpha*sss*kss^(alpha-1);
b1=1-delta+alpha*sss*kss^(alpha-1);
b2=sss*kss^(alpha-1);
b3=-css/kss;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                      Define the matrices                    %
%                                                             %
%                  F Y(t) = G Y(t+1) + H X(t+1)               %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Y :
%
% c k s
% 1 2 3
%
% X :
%
% v wc wk ws
% 1 2  3  4 
%
FF=zeros(ntots,ntots);   
GG=zeros(ntots,ntots);      
HH=zeros(ntots,ndist);

% c
FF(1,1)=-1;
GG(1,1)=-1;
GG(1,2)=a1;
GG(1,3)=a2;
HH(1,2)=-1;
HH(1,3)=a1;
HH(1,4)=a2;

% k
FF(2,1)=b3;
FF(2,2)=b1;
FF(2,3)=b2;
GG(2,2)=1;

% s
FF(3,3)=rho;
GG(3,3)=1;
HH(3,1)=-1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                        Solve the model                      %
%                                                             %
%                  F Y(t) = G Y(t+1) + H X(t+1)               %
%          ==>      Y(t) = A Y(t+1) + B X(t+1)                %
%          ==>     J Y~(t) = K Y~(t-1) + L v(t)               %
%          ==>     Y~(t) = Psi Y~(t-1) + Theta v(t)           %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AA=zeros(3,3);

AA=(FF)\(GG);

[Q,L]=eig(AA);

% Check Regularity Condition and Rearrange Q^(-1)

AMU=diag(abs(L))';
disp('Eigenvalues / modulus of the system :')
disp('=====================================')
disp(' ');       
disp([diag(L)';AMU])

[MU,k] = sort(AMU);

flag=sum(AMU<1);
if flag+npres==ntots;
   disp('==> The model is regular!')
else
   disp('==> Look out: The model is NOT regular!')
   break
end;
disp(' ');
disp('Hit a Key');pause;
disp(' ');

Q=Q(:,k);
Qs=inv(Q);

% Elements of Q^(-1) for state variables which are associated with lambda < 1

q1=Qs(1,1);
q2=Qs(1,2);
q3=Qs(1,3);

% Compute Pi

pi1=-q2/q1;
pi2=-q3/q1;

% Compute Psi and Psi

Psi  =[   b1+b3*pi1    b2+b3*pi2
          0            rho        ];

      Theta=[   0            1          ]';

disp('Policy Functions :')
disp('==================');
disp(' ');
disp('(k(t),s(t)) as functions of (k(t-1),s(t-1)) :');
disp('  ');
disp(Psi);
disp('c(t) as a function of (k(t),s(t)):');
disp(' ');
disp([pi1 pi2]);
disp('Hit a Key');pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                    IRF of k(t) and s(t)                     %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CA=[1];								    % initial value of the impulse
REPSA=zeros(npres,nrep);     				% IRF matrix of state variables
REPSA(:,1)=Theta*CA;				    % Initialize the IRF matrix of state variables
for i=2:nrep;							%
    REPSA(:,i)=Psi*REPSA(:,i-1);		% Main loop for IRF
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                        IRF of c(t)                          %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
REPXA=[  pi1   pi2   ]*REPSA;			    % IRF matrix of c(t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                            Graph IRF                        %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=1:nrep;
subplot(221);plot(T,REPXA);title('IRF(C,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(222);plot(T,REPSA(1,:));title('IRF(K,A)');
xlabel('Quarters');ylabel('% deviation')
print -dps rbcFa.eps;			     	
pause;close		                       % wait for you to hit a key and close the graph window