%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                The baseline RBC model                       %
%                                                             %
%                                                             %
%             U(C,1-h)=log(C)+theta log(1-h)                  %
%                                                             %
%                    Y=a(K^alpha)(h^(1-a))                    %
%                                                             %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;                    % Clear the memory
delete rbc.res;               % If it exists, delete the file exo.res
diary rbc.res;                % Open a new file called exo.res that will store 
										% all what is written on the screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                     Algorithm Parameters                    %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ncont=5;                  % # 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=5000;                 % # of simulations
nrep=20;                  % 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:5];             % Selection of variables of interest 
indy=3;                   % Index of Y(t) in the selection
indn=2;                   % Index of h(t) in the selection
indp=5;                   % Index of Y/N(t) in the selection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%             Structural Parameters of the Economy            %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The model is calibrated for the US econmy on Quarterly data
% The numbers are taken from Cooley and Prescott (1995) Chapter 1
%
whoy=0.6;					% Share of wages in total value added 
ykoy=3.32;					% Yearly K/Y
yiok=0.076;					% Yearly I/K
gzy=0.0156;					% Yearly per capita technological rate of growth
gny=0.0120;					% Yearly population rate of growth
hss=0.31;					% People devote 31% of their total time endowment to work
rhoa=0.95;					% Persistence of technological shock
stda=0.0079;				% Standard deviation of technological shock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                       Steady state                          %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gamy=(1+gzy)*(1+gny);						% Quarterly total rate of growth
gam=gamy^(1/4);								% Quarterly total rate of growth
alpha=1-whoy;									% elasticity of output w.r.  to capital
dy=1+yiok-(1+gzy)*(1+gny);					% Yearly depreciation rate
by=gamy/(alpha*(1/ykoy)+1-dy);			% Yearly discount factor
delta=1-(1-dy)^(1/4);						% Quarterly depreciation rate
beta=by^(1/4);									% Quarterly discount factor
yok=(gam-beta*(1-delta))/(alpha*beta);	% Model based Y/K (quarterly)
iok=gam+delta-1;								% Model based I/K (quarterly)
ioy=iok/yok;									% Model based I/Y (quarterly)
coy=1-ioy;										% Model based C/Y (quarterly)
theta=(1-hss)*(1-alpha)/(hss*coy);		% Parameter of the disutility of leisure
yss=(yok^(alpha/(alpha-1)))*hss;			% Level of output
css=coy*yss;									% Level of consumption
iss=ioy*yss;									% Level of investment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                    Matrices Coefficients                    %
%                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Static equations (c h y i p)
%
%
% Consumption
%
cc=-1;
cl=1;
%
% Hours
%
hh=1/(1-hss);
hy=-1;
hc=1;
%
% Production function
%
yy=1;
yh=alpha-1;
yk=alpha;
ya=1;
%
% Investment     
%
ic=-coy;
iy=1;
ii=-ioy;
%
% Productivity
%
py=-1;
ph=1;
pp=1;
%
% Dynamic equations (k a lambda)
%
%
% Accumulation
%
kk=gam;
kkl=delta-1;
kil=gam+delta-1;
%
% Shock
%
aa=1;
aal=-rhoa;
ae=1;      
%
% Euler
%
ll=-1;
lll=1;
lk=alpha*beta*yok/gam;
ly=lk;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%                      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 h y i p
% 1 2 3 4 5
%
% Dynamic equations :
%
% k a lambda
% 1 2 3 
%
% Mcc
%
Mcc(1,1)=cc;
Mcc(2,1)=hc;
Mcc(2,2)=hh;
Mcc(2,3)=hy;
Mcc(3,2)=yh;
Mcc(3,3)=yy;
Mcc(4,1)=ic;
Mcc(4,3)=iy;
Mcc(4,4)=ii;
Mcc(5,2)=ph;
Mcc(5,3)=py;
Mcc(5,5)=pp;
%
% Mcs
%
Mcs(1,3)=cl;
Mcs(3,1)=yk;
Mcs(3,2)=ya;
%
% Mss0
%
Mss0(1,1)=kk;
Mss0(2,2)=aa;
Mss0(3,1)=lk;
Mss0(3,3)=ll;
%
% Mss1
%
Mss1(1,1)=kkl;
Mss1(2,2)=aal;
Mss1(3,3)=lll;
%
% Msc0
%
Msc0(3,3)=ly;
%
% Msc1
%
Msc1(1,4)=kil;
%
% Mse
%
Mse(2,1)=ae;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                             %
%              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) - h(t) - Y(t) - I(t) - Y(t)/h(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(3,:));title('IRF(Y,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(223);plot(T,REPXA(4,:));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
subplot(221);plot(T,REPXA(2,:));title('IRF(H,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(222);plot(T,REPXA(5,:));title('IRF(P,A)');
xlabel('Quarters');ylabel('% deviation')
print -dps irbc2.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
%
trend=(1:long)*log(gam);
determ=[1;0;1;1;1]*trend+log([css;hss;yss;iss;yss/hss]*ones(size(trend)));
%
% stochastic Simulation 
%                   
SSIM=zeros(ntotb,slong);				% matrix of simulated state variables
SHP=zeros(nsim,length(select));		% matrix to store simulated standard deviations
SRHP=zeros(nsim,length(select));		% matrix to store simulated relative standard deviations
CHP=zeros(nsim,length(select));		% matrix to store simulated correlations
CHYN=zeros(nsim,1);						% matrix to store simulated employment-productivity correlation
%
% Beginning of Monte-Carlo loop
%
for j=1:nsim;
   disp(sprintf('Simulation : %4.0f',j));
   SHOCK=[stda]*randn(nshoc,slong);                % Draw a random number from N(0,stda)
   SSIM=MSE*SHOCK(1,1);										% Initialize each simulation
   for t=2:slong;												%
      SSIM(:,t)=MSS*SSIM(:,t-1)+MSE*SHOCK(:,t);		% Simulate state variables
   end;															%
   %
   % Retrieve control variables. Note that :
   %		1) we add the deterministic component
   %		2) we throw out the first trunc simulations in order 
   %			to be free of any initial conditions effects
   %
   XSIM=PI(select,:)*SSIM(:,trunc:trunc+long-1)+determ;
   %
   % HP-filter the control variables
	%
   HPSIM=XSIM'-HP\XSIM';
   %
   % Compute and store the moments of interest
   %
   SHP(j,:)=std(HPSIM);
   SRHP(j,:)=SHP(j,:)/std(HPSIM(:,indy));
   VCOV=cov(HPSIM);
   CHP(j,:)=VCOV(indy,:)./sqrt(diag(VCOV)'*VCOV(indy,indy));
   CHYN(j,:)=VCOV(indn,indp)./sqrt(VCOV(indn,indn)*VCOV(indp,indp));
end;                                                         
%
% End of Monte-Carlo loop
%

%
% Now, display the results
%
diary rbc.res;		% restart recording output in the file rbc.res
disp(' ')
disp('Quantitative Evaluation :');
disp('=========================');
disp(' ');
disp(sprintf('Based on %5.0f simulations',nsim));
disp(' ');
disp('Standard Deviation : C - h - Y - I - Y/N');
disp(' ');
disp(100*[mean(SHP);std(SHP)]);
disp('Relative Standard Deviation : C - h - Y - I - Y/N');
disp(' ');
disp([mean(SRHP);std(SRHP)]);
disp('Correlation with Output : C - h - Y - I - Y/N');
disp(' ');
disp([mean(CHP);std(CHP)]);
disp('Correlation Hours - productivity');
disp(' ');
disp([mean(CHYN);std(CHYN)]);
diary off;
