%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%            RBC model with fixed labor supply                %
%                                                             %
%                                                             %
%         For a description see Fabrice Collard's             %
%                manuscript with the title                    %
%  "Solving the Solow Growth Model by Linear Approximation".  %
%                                                             %
%                 This code is based on                       %
%     Fabrice's code that you can find in the appendix to     %
%  "Solving the Solow Growth Model by Linear Approximation"   %
%             and in the more recent manuscript               %
%    on the baseline RBC model (also by Fabrice Collard).     %
%                                                             %
%                Manuel Waelti, February 2002                 %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                % Clear the memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                     Algorithm Parameters                    %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ncont=3;                  % # of static equations
nback=1;                  % # of Backward endogenous state variables
nshoc=1;                  % # of exogenous shocks
ntotb=nback+nshoc;        % Total # of Backward state variables
nforw=1;                  % # of Forward endogenous state variables
nstat=ntotb+nforw;        % Total # of state variables
long=120;                 % Length of simulated series
trunc=50;                 % Truncature of simulations
slong=long+trunc;         % Length of simulation
nsim=100;                % # of simulations [original: 5000]
nrep=200;                 % IRF horizon
Mcc=zeros(ncont,ncont);   % Matrix Control-Control
Mcs=zeros(ncont,nstat);   % Matrix Control-State   
Mss0=zeros(nstat,nstat);  % Matrix State-State
Mss1=zeros(nstat,nstat);  % Matrix State-State lagged
Msc0=zeros(nstat,ncont);  % Matrix State-Control     
Msc1=zeros(nstat,ncont);  % Matrix State-Controle lagged
Mse=zeros(nstat,nshoc);   % Matrix State-Shocks           
select=[1:3];             % Selection of variables of interest [from C]
indy=2;                   % Index of Y(t) in the selection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%             Structural Parameters of the Economy            %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta=0.99;                % Discount factor
sigma=2;                  % Elasticity of utility
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                          %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
koy=(alpha*beta)/(1-beta*(1-delta));        % ratios
ioy=delta*koy;
coy=1-ioy;
yss=koy^(alpha/(1-alpha));                  % steady state values
kss=koy*yss;
css=coy*yss;
iss=ioy*yss;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                      Define the matrices                    %
%                                                             %
%                                                             %
%  Mcc X(t) = Mcs S(t)                                        %
%  Mss0 S(t+1) + Mss1 S(t) = Msc0 X(t+1) + Msc1 X(t) + e(t+1) %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Static equations :
%
% c y i
% 1 2 3
%
% Dynamic equations :
%
% k a lambda
% 1 2 3 
%
% Mcc
%
Mcc(1,1)=-sigma;
Mcc(2,2)=1;
Mcc(3,1)=-coy;
Mcc(3,2)=1;
Mcc(3,3)=-ioy;
%
% Mcs
%
Mcs(1,3)=1;
Mcs(2,1)=alpha;
Mcs(2,2)=1;
%
% Mss0
%
Mss0(1,1)=1;
Mss0(2,2)=1;
Mss0(3,1)=alpha*beta/koy;
Mss0(3,3)=-1;
%
% Mss1
%
Mss1(1,1)=-(1-delta);
Mss1(2,2)=-rho;
Mss1(3,3)=1;
%
% Msc0
%
Msc0(3,2)=alpha*beta/koy;
%
% Msc1
%
Msc1(1,3)=delta;
%
% Mse
%
Mse(2,1)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%              Solve the Farmer's system                      %
%                                                             %
%           -1                                                %
% X(t) = Mcc  Mcs S(t)                                        %
%                                                             %
%                      -1             -1                      %
% S(t+1)=(Mss0-Msc0 Mcc  Mcs)(Msc1 Mcc  Mcs-Mss1)S(t)         %
%                       -1                                    %
%        +(Mss0-Msc0 Mcc  Mcs)e(t+1)                          %
%                                                             %
% S(t+1) = W S(t) + R e(t+1)                                  %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First form the preceeding system
%
M0=inv(Mss0-Msc0*inv(Mcc)*Mcs);
M1=(Mss1-Msc1*inv(Mcc)*Mcs);
W=-M0*M1;
%
% 1) Compute the eigenvalues (MU) and eigenvectors (P) of W
% 2) Compute their modulus (AMU)
% 3) Sort the eigenvalues and rearrange the matrix of eigenvectors accordingly (P)
% 4) Form the inverse of P (Ps)
%
[P,MU] = eig(W);
AMU=diag(abs(MU))';
disp('Eigenvalues / modulus of the system :')
disp('=====================================')
disp(' ');       
disp([diag(MU)';AMU])
flag=sum(AMU>1);
if flag==nforw;
   disp('Blanchard-Kahn conditions checked')
else
   disp('Blanchard-Kahn conditions NOT checked')
end;   
disp(' ');       
[MU,k] = sort(AMU);
P=P(:,k);
Ps=inv(P);
%
% Solves the system using Farmer's method
%
M=[eye(ntotb);-inv(Ps(ntotb+1:nstat,ntotb+1:nstat))*Ps(ntotb+1:nstat,1:ntotb)];
MSS=W(1:ntotb,:)*M;
M2=M0*Mse;
MSE=M2(1:ntotb,nshoc);                                               
PI=inv(Mcc)*Mcs*M;
%
% We are done !!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%  So we now get :                                            %
%                                                             %
%  X(t)=Pi S(t)                                               %
%                                                             %
%  S(t+1)=MSS S(t) + MSE e(t+1)                               %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Policy Functions :')
disp('==================');
disp(' ');
disp('K(t+1) - A(t+1) as functions of (K(t),A(t)) :');
disp(' ');
disp(MSS);
disp('C(t) - Y(t) - I(t) as functions of (K(t),A(t)) :');
disp(' ');
disp(PI);
diary off;						% Stop recording output in a file
disp('Hit a Key');pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%     Impulse response functions to a technological shock     %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CA=[1];								    % initial value of the impulse
REPSA=zeros(ntotb,nrep);				% IRF matrix of state variables
REPSA(:,1)=MSE*CA;						% Initialize the IRF matrix of state variables
for i=2:nrep;							%
    REPSA(:,i)=MSS*REPSA(:,i-1);		% Main loop for IRF
end;									%
REPXA=PI*REPSA;							% IRF matrix of controls
%
% Graphs
%
T=1:nrep;
subplot(221);plot(T,REPXA(1,:));title('IRF(C,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(222);plot(T,REPXA(2,:));title('IRF(Y,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(223);plot(T,REPXA(3,:));title('IRF(I,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(224);plot(T,REPSA(1,:));title('IRF(K,A)');
xlabel('Quarters');ylabel('% deviation')
print -dps irbc1.eps;
pause;close		                        % wait for you to hit a key and close the graph window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                          Simulation                         %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

HP=hpmat(long,1600); % Retrieve the matrix for the HP filter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                     Loop of Simulation                      %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The deterministic component
%
% Trick : kron performs the kronecker product of 2 matrices.
%
determ=kron(log([css;yss;iss]),ones(1,long));
%
% stochastic Simulation
%
SHP=[];SRHP=[];CHP=[];
for j=1:nsim;
   disp(sprintf('Simulation : %4.0f',j));
   %
   % Random generator
   %
   CHOC=[stda]*randn(nshoc,slong);
   %
   % Initialize the states
   %
   SSIM=MSE*CHOC(1,1);
   %
   % Iterate on the state equation to obtain series of the state variables
   %
   for i=2:slong;
      SSIM=[SSIM MSS*SSIM(:,i-1)+MSE*CHOC(:,i)];
   end;
   %
   % Build selected controls
   %
   XSIM=PI(select,:)*SSIM(:,trunc:trunc+long-1)+determ;
   %
   % Obtain their HP-痩tered representation
   %
   HPSIM=XSIM'-HP\XSIM';
   %
   % Centering (Actually HP-痩tered variables are centered)
   %
   HPSIM=HPSIM-kron(mean(HPSIM),ones(length(HPSIM),1));
   %
   % Store the moments in different matrices
   %
   % For example : at iteration "i", SHP has i rows. then at iteration "i+1"
   % the SHP equals the SHP matrix to which we add another
   % line containing std(HPSIM).
   %
   SHP=[SHP;std(HPSIM)];
   SRHP=[SRHP;std(HPSIM)/std(HPSIM(:,indy))];
   VCOV=HPSIM'*HPSIM;
   CHP=[CHP;VCOV(indy,:)./sqrt(diag(VCOV)'*VCOV(indy,indy))];
end;                                                         
%
% End of Monte-Carlo loop. Now, display the results
%
diary solow_new.res;		% restart recording output in the file solow_new.res
disp(' ');
disp('Quantitative Evaluation :');
disp('=========================');
disp(' ');
disp(sprintf('%Based on %5.0f simulations',nsim));
disp(' ');
disp('Standard Deviations : C - Y - I');
disp(' ');
disp(100*[mean(SHP);std(SHP)]);
disp(' ');
disp('Relative Standard Deviations : C - Y - I');
disp(' ');
disp([mean(SRHP);std(SRHP)]);
disp(' ');
disp('Correlation with Output : C - Y - I');
disp(' ');
disp([mean(CHP);std(CHP)]);
diary off;