%% Phase Diagram Extended Ramsey growth model

% Last edited by Ben on 29 September 2015

% Endogenous labour supply (unit elastic case)
% Parameters as chosen in equation (13.119) in FoMM3

clear
clf
clc
close all

% Structural parameters
PAR.alpha = 1/3;              % Efficiency parameter capital
PAR.eps_C = 0.183;            % Utility parameter of consumption
PAR.delta = (1.1)^(1/4)-1;    % Quarterly depreciation rate for capital
PAR.rho   = (1.065)^(1/4)-1;  % Quarterly rate of time preference


% Exogenous variables
EXO.G     = 0.20;             % Government consumption
%EXO.G     = 0;
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        = PAR.rho;
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
% Set the scale parameter production function to get Y0 right
PAR.Z_0      = K0^(-PAR.alpha) * L0^(PAR.alpha-1) * Y0;     
% Set the taste parameter eps_C to get L0 right
PAR.om_LL    = (1 - L0)/L0 ;
PAR.om_C     = C0/Y0;
PAR.om_G     = EXO.G/Y0;
PAR.eps_C    = PAR.om_C / (PAR.om_C + (1 - PAR.alpha) * PAR.om_LL);
Gamma_0      = PAR.Z_0 * ((1 - PAR.alpha * PAR.eps_C) /(1 - PAR.eps_C)) ;

% Some critical values
%PAR.L_min    = PAR.eps_C * (1-PAR.alpha) /(PAR.eps_C * (1-PAR.alpha) + (1-PAR.eps_C) * (1- PAR.om_G) );
PAR.L_min     = PAR.eps_C * (1-PAR.alpha) /(1 - PAR.alpha * PAR.eps_C) ;
PAR.K_crit    = PAR.alpha * EXO.G / ((1-PAR.alpha)*PAR.delta) ;
funmincapital = @(x) (PAR.delta * x + EXO.G) .* x.^(-PAR.alpha) - PAR.Z_0 ;
PAR.K_min     = fzero(funmincapital,0.01);
PAR.K_max     = fzero(funmincapital,500);
F_min         = (EXO.G/(1-PAR.alpha))^(1-PAR.alpha) * (PAR.delta/PAR.alpha)^PAR.alpha ;
funminlabour  = @(x) Gamma_0 * (x - PAR.L_min) .* x.^(-PAR.alpha) - F_min;          
PAR.L_crit    = fzero(funminlabour,0.2)
        
        
% Auxiliary functions f(L) and g(K,G)
L_vals    = linspace(PAR.L_min,1,1000);
K_vals    = linspace(PAR.K_min,PAR.K_max,5000);
Gfunction = Gamma_0 * (L_vals - PAR.L_min) .* L_vals.^(-PAR.alpha);  
Ffunction = (PAR.delta * K_vals + EXO.G) .* K_vals.^(-PAR.alpha);
% F_min     = min(Ffunction)
Ben = PAR.Z_0 * PAR.K_min^PAR.alpha - PAR.delta * PAR.K_min - EXO.G  ;        

% Consumption isocline
Ciso = (PAR.rho+PAR.delta) * ((1 - PAR.alpha) * PAR.eps_C) /(PAR.alpha*(1 - PAR.eps_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.5])
xlabel('capital')
ylabel('\phi(K,G) function')
plot(K_vals,Ffunction,'-','LineWidth',3)
plot([0 450],[0.5 0.5],':','LineWidth',1)
plot([0 450],[PAR.Z_0 PAR.Z_0],':','LineWidth',1)
plot([81 81],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13A_02_phi_K') 
hold off
pause(1)

figure(2)
hold on
box off
ylim([0,1.5])
xlabel('labour')
ylabel('\xi(L) function')
plot(L_vals,Gfunction,'-','LineWidth',3)
plot([0 1],[0.5 0.5],':','LineWidth',1)
plot([0 1],[PAR.Z_0 PAR.Z_0],':','LineWidth',1)
plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13A_02_xi_L') 
hold off
pause(1)


% Compute equilibrium L for each K value in [K_min,K_max]
        Lpad          = zeros(size(K_vals));
        Cpad          = zeros(size(K_vals));
        
        L0   = 0.9;
        
        for i=1:numel(K_vals)
        funeqlab  = @(x) (PAR.delta * K_vals(i) + EXO.G) .* x.^(PAR.alpha) ...
           - Gamma_0 * (x - PAR.L_min) .* K_vals(i).^(PAR.alpha);          
        Lpad(i)   = fzero(funeqlab,L0);
        Cpad(i)   = PAR.Z_0 * Lpad(i).^(1-PAR.alpha) .* K_vals(i).^(PAR.alpha) ...
            - EXO.G - PAR.delta * K_vals(i);
        L0 = Lpad(i);
        end


        
figure(3)
hold on
box off
%xlim([0,10])
ylim([0,1])
xlabel('capital stock K(t)')
ylabel('consumption C(t)')
plot(K_vals,Cpad,'-','LineWidth',3)
plot(K_vals,Ciso,'-','LineWidth',3)
text(45,0.2,'C isocline')
text(225,0.7,'K isocline')
hold off
print('-deps2','-loose','F13A_03') 
pause(3)


etagamma = (1/Gamma_0) *(Lpad./K_vals).^(PAR.alpha) .* ...
     ((1-PAR.alpha)*PAR.delta*K_vals - PAR.alpha*EXO.G) ...
      ./ ( (1-PAR.alpha)*Lpad + PAR.alpha * PAR.L_min ) ;

Kisoslope = PAR.Z_0 * (PAR.alpha + (1-PAR.alpha)* etagamma) ...
    .* (Lpad./K_vals).^(1-PAR.alpha) - PAR.delta ;
  
figure(4)
hold on
box off
%xlim([0,10])
ylim([0,1])
xlabel('capital stock K(t)')
ylabel('\eta_{\gamma}(K)')
plot(K_vals,etagamma,'-','LineWidth',3)
% text(45,0.2,'C isocline')
% text(225,0.7,'K isocline')
hold off
print('-deps2','-loose','F13A_04') 
pause(3)
      
figure(5)
hold on
box off
xlim([0,20])
%ylim([0,1])
xlabel('capital stock K(t)')
ylabel('slope K isocline')
plot(K_vals,Kisoslope,'-','LineWidth',3)
% text(45,0.2,'C isocline')
% text(225,0.7,'K isocline')
hold off
print('-deps2','-loose','F13A_05') 
pause(3)
    