%% Extended Blanchard-Yaari model
%
% Last edited by Ben on 29 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); 

% Compute the steady-state equilibrium

x = zeros(7,1);
[x(:,1),fval,exitflag,output] = fsolve(@Function15_02_SS,[W0,R0,C0,K0,Y0,I0,L0]'...
    ,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-6,'TolFun',1.0e-6),PAR,EXO);
                                     
fval;
output;
x;
VAR.Wbase  = x(1,1);
VAR.Rbase  = x(2,1);
VAR.Cbase  = x(3,1);
VAR.Kbase  = x(4,1);
VAR.Ybase  = x(5,1);
VAR.Ibase  = x(6,1);
VAR.Lbase  = x(7,1);
VAR.TRbase = EXO.t_L * VAR.Wbase * VAR.Lbase + EXO.t_C * VAR.Cbase;

% Tax shock in the macro model: t_C from 0.1 to 0.2
EXO.t_C   = 0.20 ;

x = zeros(7,1);
[x(:,1),fval,exitflag,output] = fsolve(@Function15_02_SS,[W0,R0,C0,K0,Y0,I0,L0]'...
    ,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-6,'TolFun',1.0e-6),PAR,EXO);
                                     
fval;
output;
x;
VAR.Wt_Cinc = x(1,1);
VAR.Rt_Cinc = x(2,1);
VAR.Ct_Cinc = x(3,1);
VAR.Kt_Cinc = x(4,1);
VAR.Yt_Cinc = x(5,1);
VAR.It_Cinc = x(6,1);
VAR.Lt_Cinc = x(7,1);

EXO.t_C     = 0.10 ;


% Model parameters
    disp(' ')
    disp('Model parameters')
    disp(sprintf('      Efficiency parameter capital (alpha):                   %12.8f',PAR.alpha))
    disp(sprintf('      Scale parameter production function (Z_0):              %12.8f',PAR.Z_0))
    disp(sprintf('      Quarterly capital depreciation rate (delta):            %12.8f',PAR.delta))
    disp(sprintf('      Pure rate of time preference (rho):                     %12.8f',PAR.rho))
    disp(sprintf('      Mortality rate (mu):                                    %12.8f',PAR.mu))
    disp(' ')    
    disp('Exogenous variables')
    disp(sprintf('      Labour income tax (t_L):                                %12.8f',EXO.t_L))
    disp(sprintf('      Consumption (t_C):                                      %12.8f',EXO.t_C))
    disp(sprintf('      Government consumption (G):                             %12.8f',EXO.G))

% Implied calibration
    disp(' ')
    disp('Steady-state equilibrium')
    disp(sprintf('      Consumption (C):                                        %12.8f',VAR.Cbase))
    disp(sprintf('      Capital stock (K):                                      %12.8f',VAR.Kbase))
    disp(sprintf('      Employment (L):                                         %12.8f',VAR.Lbase))
    disp(sprintf('      Output (Y):                                             %12.8f',VAR.Ybase))
    disp(sprintf('      Investment (I):                                         %12.8f',VAR.Ibase))
    disp(sprintf('      Interest rate (R):                                      %12.8f',VAR.Rbase))
    disp(sprintf('      Wage rate (W):                                          %12.8f',VAR.Wbase))
    disp(' ')
    disp('Consumption tax increase from t_C=0.1 to t_C=0.2')
    disp(sprintf('      Consumption (C):                                        %12.8f',VAR.Ct_Cinc))
    disp(sprintf('      Capital stock (K):                                      %12.8f',VAR.Kt_Cinc))
    disp(sprintf('      Employment (L):                                         %12.8f',VAR.Lt_Cinc))
    disp(sprintf('      Output (Y):                                             %12.8f',VAR.Yt_Cinc))
    disp(sprintf('      Investment (I):                                         %12.8f',VAR.It_Cinc))
    disp(sprintf('      Interest rate (R):                                      %12.8f',VAR.Rt_Cinc))
    disp(sprintf('      Wage rate (W):                                          %12.8f',VAR.Wt_Cinc))
    

% Solve the model in continuous time using BVP4C 
    
% Shock time:            t_A = 0
% Implementation time:   t_I
% End time:              t_E
% Unanticipated / permanent shock t_I = 0, t_E = inf
% Anticipated   / permanent shock t_I > 0, t_E = inf
% Unanticipated / permanent shock t_I = 0, t_E = finite
% Anticipated   / permanent shock t_I > 0, t_E = finite




    CaseVal = 2;
        
    switch CaseVal
        case 1
        % Convergence to the steady state following a twenty percent
        % destruction of the capital stock
        
        VAR.K_initial = 0.8 * VAR.Kbase;
        VAR.Keq0      = VAR.K_initial ;           
        VAR.Keq1      = VAR.Kbase ; 
        VAR.Ceq1      = VAR.Cbase ; 

        [time,Cpath,Kpath] = Function15_02_TP(PAR,EXO,VAR) ;
        
        % Compute equilibrium employment
        Lpath         = zeros(size(Kpath));
       
        
        for i=1:length(Kpath)
        funeqemp  = @(x) (1 - PAR.alpha) * (1 - EXO.t_L) * PAR.Z_0 * Kpath(i)^PAR.alpha * (1 - x)  ...
          - ((1-PAR.epsilon)/PAR.epsilon) * (1 + EXO.t_C) * Cpath(i) * x^PAR.alpha ;
        Lpath(i)   = fzero(funeqemp,0.3);
        end
        
        Ypath     = PAR.Z_0 * Kpath.^PAR.alpha .* Lpath.^(1-PAR.alpha) ;
        Wpath     = (1-PAR.alpha) * Ypath ./ Lpath ;
        Rpath     = PAR.alpha * (Ypath ./ Kpath) - PAR.delta ;
        Ipath     = Ypath - Cpath - EXO.G ;
        TRpath    = EXO.t_L * Wpath .* Lpath + EXO.t_C * Cpath;
        

figure(1)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('capital')
plot(time,Kpath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_K') 
pause(1)

figure(2)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('consumption')
plot(time,Cpath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_C') 
pause(1)

figure(3)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('employment')
plot(time,Lpath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_L') 
pause(1)

figure(4)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('output')
plot(time,Ypath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_Y') 
pause(1)

figure(5)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('transfers')
plot(time,TRpath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_TR') 
pause(1)

figure(6)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('wages')
plot(time,Wpath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_W') 
pause(1)

figure(7)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('interest rate')
plot(time,Rpath,'k','LineWidth',3);
hold off
%print('-deps2','-loose','Figure_ExtendedBlanchard__Case1_R') 
pause(1)

        case 2
        % Transition path following an unanticipated and permanent increase 
        % in the consumption tax from t_C = 0.1 to t_C = 0.2
        
        EXO.t_C       = 0.2;
        VAR.Keq0      = VAR.Kbase ;           
        VAR.Keq1      = VAR.Kt_Cinc ; 
        VAR.Ceq1      = VAR.Ct_Cinc ; 

        [time,Cpath,Kpath] = Function15_02_TP(PAR,EXO,VAR) ;
        
        % Compute equilibrium employment
        Lpath         = zeros(size(Kpath));
       
        
        for i=1:length(Kpath)
        funeqemp  = @(x) (1 - PAR.alpha) * (1 - EXO.t_L) * PAR.Z_0 * Kpath(i)^PAR.alpha * (1 - x)  ...
          - ((1-PAR.epsilon)/PAR.epsilon) * (1 + EXO.t_C) * Cpath(i) * x^PAR.alpha ;
        Lpath(i)   = fzero(funeqemp,0.3);
        end
        
        Ypath     = PAR.Z_0 * Kpath.^PAR.alpha .* Lpath.^(1-PAR.alpha) ;
        Wpath     = (1-PAR.alpha) * Ypath ./ Lpath ;
        Rpath     = PAR.alpha * (Ypath ./ Kpath) - PAR.delta ;
        Ipath     = Ypath - Cpath - EXO.G ;
        TRpath    = EXO.t_L * Wpath .* Lpath + EXO.t_C * Cpath;
        
        EXO.t_C       = 0.1;

figure(1)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('capital')
plot(time,Kpath,'k','LineWidth',3);
plot([0 120],[VAR.Kbase VAR.Kbase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_K') 
pause(1)

figure(2)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('consumption')
plot(time,Cpath,'k','LineWidth',3);
plot([0 120],[VAR.Cbase VAR.Cbase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_C') 
pause(1)

figure(3)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('employment')
plot(time,Lpath,'k','LineWidth',3);
plot([0 120],[VAR.Lbase VAR.Lbase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_L') 
pause(1)

figure(4)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('output')
plot(time,Ypath,'k','LineWidth',3);
plot([0 120],[VAR.Ybase VAR.Ybase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_Y') 
pause(1)

figure(5)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('transfers')
plot(time,TRpath,'k','LineWidth',3);
plot([0 120],[VAR.TRbase VAR.TRbase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_TR') 
pause(1)

figure(6)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('wages')
plot(time,Wpath,'k','LineWidth',3);
plot([0 120],[VAR.Wbase VAR.Wbase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_W') 
pause(1)

figure(7)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('interest rate')
plot(time,Rpath,'k','LineWidth',3);
plot([0 120],[VAR.Rbase VAR.Rbase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_R') 
pause(1)

figure(8)
hold on
box off
xlim([0,120])
xlabel('post-shock time (quarters)')
ylabel('investment')
plot(time,Ipath,'k','LineWidth',3);
plot([0 120],[VAR.Ibase VAR.Ibase],'--','LineWidth',2)
hold off
print('-deps2','-loose','F15_19_I') 
pause(1)



     end
      
    