%% Phase Diagram Extended Blanchard model

% Last edited by Ben on 28 October 2015

% Endogenous labour supply (unit elastic case)
% Parameters as chosen in Section 15.5.4 in FoMM3

clear
clf
clc
close all

% Structural parameters
PAR.alpha     = 1/3;                  % Efficiency parameter capital
PAR.delta     = (1.1)^(1/4)-1;        % Quarterly depreciation rate for capital
PAR.mu        = 0.005;                % Quarterly mortality rate

% Exogenous variables
EXO.t_L       = 0.30;                 % Labour income tax
EXO.t_C       = 0.10;                 % Consumption tax
EXO.G         = 0.00;                 % Government consumption
EXO.time      = linspace(0,200,2001); % Post-shock time (for plotting) 

% Target variables
Y0        = 1.00;                     % Base line per capita output per quarter 
L0        = 0.20;
W0        = (1-PAR.alpha) * Y0 / L0;
R0        = (1.065)^(1/4)-1; 
KYratio0  = PAR.alpha/(R0+PAR.delta); % Base line capital-output ratio
K0        = KYratio0 * Y0;
I0        = PAR.delta * K0;
C0        = Y0 - I0 - EXO.G; 

% Calibration parameters
% Technology index
PAR.Z_0      = K0^(-PAR.alpha) * L0^(PAR.alpha-1) * Y0;    
% Utility parameter of consumption
PAR.epsilon  = (1 + ( (1-PAR.alpha) * (1-L0)/L0 * (1-EXO.t_L) ) ...
    / (C0/Y0 * (1+EXO.t_C) ))^(-1);
 % Quarterly rate of time preference
gamma_0      = PAR.epsilon * K0/C0 / (1+EXO.t_C) ; 
PAR.rho      = (R0 - PAR.mu^2 * gamma_0) / (1 + PAR.mu * gamma_0); 


% Some critical values
PAR.L_min     = PAR.epsilon * (1-PAR.alpha) *(1 - EXO.t_L) ...
    /( PAR.epsilon * (1-PAR.alpha) *(1 - EXO.t_L) + (1 - PAR.epsilon)*(1+EXO.t_C)) ;
PAR.L_max     = 1;
PAR.K_min     = 0;
PAR.K_max     = (PAR.Z_0/PAR.delta)^(1/(1-PAR.alpha));
        
    

L_vals    = linspace(PAR.L_min,1,1000);
K_vals    = linspace(PAR.K_min,PAR.K_max,5000);
Gamma_0   = PAR.epsilon * (1-PAR.alpha) *(1 - EXO.t_L) ...
    / (PAR.delta * (1 - PAR.epsilon) * (1 + EXO.t_C)) ;
K_kiso = (PAR.Z_0 * (Gamma_0/PAR.L_min) * (L_vals - PAR.L_min) .* L_vals.^(-PAR.alpha)).^(1/(1-PAR.alpha)) ;  
C_kiso = PAR.Z_0 * K_kiso.^(PAR.alpha) .* L_vals.^(1-PAR.alpha) - PAR.delta * K_kiso ;
%plot(K_kiso,C_kiso)

L_space    = linspace(0.01,1,1000);
zeta_0 = PAR.mu * (PAR.rho+PAR.mu)*(1-PAR.epsilon) ...
    / (PAR.alpha *(1-PAR.alpha) * (1+EXO.t_C) ) ;
ystar = (PAR.rho+PAR.delta) / PAR.alpha ;

y = ystar * (1/2 + sqrt(1/4 + zeta_0 * ystar^(-2) * L_space ./ (1-L_space) ));
KLratio = (PAR.Z_0 ./ y).^(1/(1-PAR.alpha)) ;

K_ciso = KLratio .* L_space;
C_ciso = (PAR.epsilon * (1-PAR.alpha) *(1 - EXO.t_L) * PAR.Z_0 ...
    / ((1 - PAR.epsilon) * (1 + EXO.t_C))) * KLratio.^PAR.alpha .* (1 - L_space);
plot(K_ciso,C_ciso)

% Consumption isocline RA model
Ciso = ((PAR.rho+PAR.delta) * (1-EXO.t_L) * (1 - PAR.alpha) * PAR.epsilon) ...
    /(PAR.alpha*(1 - PAR.epsilon)*(1 + EXO.t_C)) ...
    * ( (PAR.alpha * PAR.Z_0 /(PAR.rho+PAR.delta))^(1/(1-PAR.alpha)) - K_vals);

figure(1)
hold on
box off
ylim([0,1.2])
xlim([0,50])
xlabel('capital stock K(t)')
ylabel('consumption C(t)')
plot(K_kiso,C_kiso,'-','LineWidth',3)
plot(K_ciso,C_ciso,'-','LineWidth',3)
plot(K_vals,Ciso,'--','LineWidth',2)
text(20,0.17,'C isocline','FontSize',12)
text(35,1.00,'K isocline','FontSize',12)
print('-deps2','-loose','F15A_01') 
hold off
pause(1)

    