% This Matlab m file actually plot the phase diagram for classical RBC model.
% This code is written by Xiangyang Li, Ph.D in SUFE, School of Finance
% xli21@nd.edu
clear;
%set in quarterly frequency
alpha=.36;
beta =.99;
sigma =2;
delta=.025;
%steady states and the golden capital stock.
ks = (alpha/(1/beta - 1 + delta))^(1/(1-alpha));
kgold = (alpha/delta)^(1/(1-alpha));
cs = ks^alpha - delta*ks;
%draw the isocline
kv=0:.1:200;
cv=0:.1:6;
% the capital isocline
plot(kv,kv.^alpha-delta*kv,'-r','LineWidth',2);
hold on;
% the consumption isocline
plot(ks*ones(size(cv)),cv,'-b','LineWidth',2);
plot(kgold*ones(size(cv)),cv,'-g');
%1. suppose we start from the point c_t = 2 and k_t =20, to see what happend
T = 52;
C=zeros(1,T);
K=zeros(1,T);
C(1)=2;
K(1) = 20;
for ii = 1:T
    %resources constraint
    K(ii+1) = K(ii)^alpha + (1 -delta)*K(ii)  - C(ii);
    %FOC of the problem
    C(ii+1) =( C(ii)^(-sigma)/(alpha*K(ii+1)^(alpha -1) + (1-delta) ))^(-1/sigma);
end
plot(K,C,'-k');
%2. suppose we start from the point c_t = 1 and k_t =20, to see what happend
T = 120;
C=zeros(1,T);
K=zeros(1,T);
C(1)=1;
K(1) = 20;
for ii = 1:T
    K(ii+1) = K(ii)^alpha + (1 -delta)*K(ii)  - C(ii);
    C(ii+1) =( C(ii)^(-sigma)/(alpha*K(ii+1)^(alpha -1) + (1-delta) ))^(-1/sigma);
end
plot(K,C,'-m');
%3. suppose we start from the point c_t = 4.5 and k_t =70, to see what happend
T = 70;
C=zeros(1,T);
K=zeros(1,T);
C(1)=4.5;
K(1) = 120;
for ii = 1:T
    K(ii+1) = K(ii)^alpha + (1 -delta)*K(ii)  - C(ii);
    C(ii+1) =( C(ii)^(-sigma)/(alpha*K(ii+1)^(alpha -1) + (1-delta) ))^(-1/sigma);
end
plot(K,C,'-c');
%4. suppose we start from the point c_t = 3.5 and k_t =90, to see what happend
T = 370;
C=zeros(1,T);
K=zeros(1,T);
C(1)=3.1;
K(1) = 80;
for ii = 1:T
    K(ii+1) = K(ii)^alpha + (1 -delta)*K(ii)  - C(ii);
    C(ii+1) =( C(ii)^(-sigma)/(alpha*K(ii+1)^(alpha -1) + (1-delta) ))^(-1/sigma);
end
plot(K,C,'-r');

%5. plot the policy function, policy function coefficient is obtained
%by solving the RBC model.
pol=0.462886778502290;
K=10:0.1:110;
C=pol*K*cs/ks + cs*(1-pol);
plot(K,C,'-k','LineWidth',2);
hold off;
legend('captial isocline', 'consumption isocline', ...
    'capital gold','explosive path 1', 'explosive path 2',...
'explosive path 3', 'explosive path 4','Saddle Path');