%% Unit Elastic RCK Model
%
%  Figure 13.8: Permanent or temporary fiscal policy
%
% Ben J. Heijdra, Foundations of Modern Macroeconomics, Third Ed.
% Oxford University Press
%
% Chapter 13
%
% Ben J. Heijdra
% Groningen, September 24, 2015



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
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.eps_C    = PAR.om_C / (PAR.om_C + (1-PAR.alpha) * PAR.om_LL);
PAR.phi      = (1 + PAR.om_LL) / (1 + PAR.alpha * PAR.om_LL);

Y0 = 0.95;

% Compute the initial steady-state equilibrium

x = zeros(7,1);
[x(:,1),fval,exitflag,output] = fsolve(@StStMacroModelLabSup,[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.Wstar = x(1,1);
VAR.Rstar = x(2,1);
VAR.Cstar = x(3,1);
VAR.Kstar = x(4,1);
VAR.Ystar = x(5,1);
VAR.Istar = x(6,1);
VAR.Lstar = x(7,1);

% Compute the new steady-state equilibrium
% G(t) from G_0=0.20 to G_1=0.25

dG        = 0.05 ;
EXO.G     = EXO.G + dG;                 % New government consumption
x = zeros(7,1);
[x(:,1),fval,exitflag,output] = fsolve(@StStMacroModelLabSup,[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.Winf = x(1,1);
VAR.Rinf = x(2,1);
VAR.Cinf = x(3,1);
VAR.Kinf = x(4,1);
VAR.Yinf = x(5,1);
VAR.Iinf = x(6,1);
VAR.Linf = x(7,1);

EXO.G     = EXO.G - dG;                 % Reset government consumption

%  Exact long-run effects
dYdGinf   = (VAR.Yinf - VAR.Ystar) / dG 
dCdGinf   = (VAR.Cinf - VAR.Cstar) / dG 
dIdGinf   = (VAR.Iinf - VAR.Istar) / dG 
ElastKG   = ((VAR.Kinf - VAR.Kstar)/VAR.Kstar) / (dG/EXO.G)
ElastLG   = ((VAR.Linf - VAR.Lstar)/VAR.Lstar) / (dG/EXO.G)


%
%  Characteristic roots of the Jacobian matrix
%
delta_11   = -(1-PAR.alpha) * PAR.om_LL * (PAR.rho + PAR.delta) ...
    / (1 + PAR.alpha * PAR.om_LL) 
delta_12   = -(1-PAR.alpha) * (PAR.rho + PAR.delta) * (VAR.Cstar/VAR.Kstar) ...
    / (1 + PAR.alpha * PAR.om_LL) 
delta_21   = -1 -  (1-PAR.alpha) * PAR.om_LL * (VAR.Ystar/VAR.Cstar) ...
    / (1 + PAR.alpha * PAR.om_LL) 
delta_22   = (VAR.Cstar+EXO.G)/VAR.Kstar - (1-PAR.alpha) * (VAR.Ystar/VAR.Kstar) ...
    / (1 + PAR.alpha * PAR.om_LL) 
%
%  Determinant and characteristic roots
%
DDelta    =  delta_11*delta_22 - delta_12*delta_21
TDelta    =  delta_11 + delta_22 
XX3       =  sqrt(TDelta*TDelta - 4*DDelta) ;
lambda_1  =  (-TDelta + XX3) / 2;
lambda_2  =  ( TDelta + XX3) / 2;


% Plot figures for different values of xi_G 

xi_G1  = 0;
xi_G2  = 0.02;
xi_G3  = 0.10;

Gpath1 = exp(-xi_G1 * time)  * dG;
Gpath2 = exp(-xi_G2 * time)  * dG;
Gpath3 = exp(-xi_G3 * time)  * dG;

figure(1)
hold on
box off
ylim([0,0.05])
% xlabel('postshock time (quarters)')
% ylabel('G(t)-G_0')
plot(time,Gpath1,'-','LineWidth',3)
plot(time,Gpath2,'--','LineWidth',3)
plot(time,Gpath3,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
 print('-deps2','-loose','F13_08_Gpath') 
hold off
pause(1)


Transpath1 = (exp(-lambda_1 * time) - exp(-xi_G1 * time))  ...
    / (xi_G1 - lambda_1) ;
Transpath2 = (exp(-lambda_1 * time) - exp(-xi_G2 * time))  ...
    / (xi_G2 - lambda_1);
Transpath3 = (exp(-lambda_1 * time) - exp(-xi_G3 * time))  ...
    / (xi_G3 - lambda_1);

Tdotpath1 = (lambda_1 * exp(-lambda_1 * time) - xi_G1 * exp(-xi_G1 * time)) ...
    / (lambda_1 - xi_G1) ;
Tdotpath2 = (lambda_1 * exp(-lambda_1 * time) - xi_G1 * exp(-xi_G2 * time)) ...
    / (lambda_1 - xi_G2) ;
Tdotpath3 = (lambda_1 * exp(-lambda_1 * time) - xi_G1 * exp(-xi_G3 * time)) ...
    / (lambda_1 - xi_G3) ;

figure(2)
hold on
box off
%ylim([0,0.1])
% xlabel('postshock time (quarters)')
% ylabel('T(\xi_G,\lambda_1,t) dG')
plot(time,lambda_1 * Transpath1,'-','LineWidth',3)
plot(time,lambda_1 * Transpath2,'--','LineWidth',3)
plot(time,lambda_1 * Transpath3,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Transpath') 
hold off
pause(1)

return 

Czero1 = VAR.Cstar + (lambda_2 - delta_11) * dG ...
    /(delta_21 * (lambda_2 + xi_G1)) ;
Czero2 = VAR.Cstar + (lambda_2 - delta_11) * dG ...
    /(delta_21 * (lambda_2 + xi_G2)) ;
Czero3 = VAR.Cstar + (lambda_2 - delta_11) * dG ...
    /(delta_21 * (lambda_2 + xi_G3)) ;

Cpath1 = Czero1 * exp(-lambda_1 .* time) ...
    + VAR.Cstar * (1 - exp(-lambda_1 .* time)) ...
    + delta_12 /(lambda_2 + xi_G1) * dG * Transpath1 ;
Cpath2 = Czero2 * exp(-lambda_1 .* time) ...
    + VAR.Cstar * (1 - exp(-lambda_1 .* time)) ...
    + delta_12 /(lambda_2 + xi_G2) * dG * Transpath2 ;
Cpath3 = Czero3 * exp(-lambda_1 .* time) ...
    + VAR.Cstar * (1 - exp(-lambda_1 .* time)) ...
    + delta_12 /(lambda_2 + xi_G3) * dG * Transpath3 ;


Cdotpath1 = lambda_1 * (VAR.Cstar - Czero1) * exp(-lambda_1 .* time) ...
     + delta_12 /(lambda_2 + xi_G1) * dG *  Tdotpath1 ;
Cdotpath2 = lambda_1 * (VAR.Cstar - Czero1) * exp(-lambda_1 .* time) ...
     + delta_12 /(lambda_2 + xi_G2) * dG *  Tdotpath2 ;
Cdotpath3 = lambda_1 * (VAR.Cstar - Czero1) * exp(-lambda_1 .* time) ...
     + delta_12 /(lambda_2 + xi_G3) * dG *  Tdotpath3 ;


figure(3)
hold on
box off
%ylim([0,0.1])
% xlabel('postshock time (quarters)')
% ylabel('(C(t)-C^*)/dG')
plot(time,(Cpath1-VAR.Cstar)/dG,'-','LineWidth',3)
plot(time,(Cpath2-VAR.Cstar)/dG,'--','LineWidth',3)
plot(time,(Cpath3-VAR.Cstar)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Cpath') 
hold off
pause(1)


Kpath1 = VAR.Kstar - (delta_11 + xi_G1) /(lambda_2 + xi_G1) * dG * Transpath1 ;
Kpath2 = VAR.Kstar - (delta_11 + xi_G2) /(lambda_2 + xi_G2) * dG * Transpath2 ;
Kpath3 = VAR.Kstar - (delta_11 + xi_G3) /(lambda_2 + xi_G3) * dG * Transpath3 ;
Kdotpath1 = - (delta_11 + xi_G1) /(lambda_2 + xi_G1) * dG * Tdotpath1 ;
Kdotpath2 = - (delta_11 + xi_G2) /(lambda_2 + xi_G2) * dG * Tdotpath2 ;
Kdotpath3 = - (delta_11 + xi_G3) /(lambda_2 + xi_G3) * dG * Tdotpath3 ;


figure(4)
hold on
box off
%ylim([0,0.1])
% xlabel('postshock time (quarters)')
% ylabel('(K(t)-K^*)/dG')
plot(time,(EXO.G/VAR.Kstar)*(Kpath1-VAR.Kstar)/dG,'-','LineWidth',3)
plot(time,(EXO.G/VAR.Kstar)*(Kpath2-VAR.Kstar)/dG,'--','LineWidth',3)
plot(time,(EXO.G/VAR.Kstar)*(Kpath3-VAR.Kstar)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Kpath') 
hold off
pause(1)


Ypath1 = VAR.Ystar + PAR.alpha * PAR.phi * (VAR.Ystar/VAR.Kstar) * ...
    (Kpath1 - VAR.Kstar) - (PAR.phi - 1) * (VAR.Ystar/VAR.Cstar) * ...
    (Cpath1 - VAR.Cstar) ;
Ypath2 = VAR.Ystar + PAR.alpha * PAR.phi * (VAR.Ystar/VAR.Kstar) * ...
    (Kpath2 - VAR.Kstar) - (PAR.phi - 1) * (VAR.Ystar/VAR.Cstar) * ...
    (Cpath2 - VAR.Cstar) ;
Ypath3 = VAR.Ystar + PAR.alpha * PAR.phi * (VAR.Ystar/VAR.Kstar) * ...
    (Kpath3 - VAR.Kstar) - (PAR.phi - 1) * (VAR.Ystar/VAR.Cstar) * ...
    (Cpath3 - VAR.Cstar) ;
Ydotpath1 = PAR.alpha * PAR.phi * (VAR.Ystar/VAR.Kstar) * Kdotpath1 ...
     - (PAR.phi - 1) * (VAR.Ystar/VAR.Cstar) * Cdotpath1 ;
Ydotpath2 = PAR.alpha * PAR.phi * (VAR.Ystar/VAR.Kstar) * Kdotpath2 ...
     - (PAR.phi - 1) * (VAR.Ystar/VAR.Cstar) * Cdotpath2 ;
Ydotpath3 = PAR.alpha * PAR.phi * (VAR.Ystar/VAR.Kstar) * Kdotpath3 ...
     - (PAR.phi - 1) * (VAR.Ystar/VAR.Cstar) * Cdotpath3 ;


figure(5)
hold on
box off
ylim([0,1.2])
% xlabel('postshock time (quarters)')
% ylabel('(Y(t)-Y^*)/dG')
plot(time,(Ypath1-VAR.Ystar)/dG,'-','LineWidth',3)
plot(time,(Ypath2-VAR.Ystar)/dG,'--','LineWidth',3)
plot(time,(Ypath3-VAR.Ystar)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Ypath') 
hold off
pause(1)


Lpath1 = VAR.Lstar + (PAR.phi - 1) * (VAR.Lstar/(1-PAR.alpha)) * ...
    ( PAR.alpha * (Kpath1 - VAR.Kstar)/VAR.Kstar ...
                - (Cpath1 - VAR.Cstar)/VAR.Cstar  );
Lpath2 = VAR.Lstar + (PAR.phi - 1) * (VAR.Lstar/(1-PAR.alpha)) * ...
    ( PAR.alpha * (Kpath2 - VAR.Kstar)/VAR.Kstar ...
                - (Cpath2 - VAR.Cstar)/VAR.Cstar  );
Lpath3 = VAR.Lstar + (PAR.phi - 1) * (VAR.Lstar/(1-PAR.alpha)) * ...
    ( PAR.alpha * (Kpath3 - VAR.Kstar)/VAR.Kstar ...
                - (Cpath3 - VAR.Cstar)/VAR.Cstar  );
            
figure(6)
hold on
box off
%ylim([0,1.2])
% xlabel('postshock time (quarters)')
% ylabel('(L(t)-L^*)/dG')
plot(time,(EXO.G/VAR.Lstar)*(Lpath1-VAR.Lstar)/dG,'-','LineWidth',3)
plot(time,(EXO.G/VAR.Lstar)*(Lpath2-VAR.Lstar)/dG,'--','LineWidth',3)
plot(time,(EXO.G/VAR.Lstar)*(Lpath3-VAR.Lstar)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Lpath') 
hold off
pause(1)

Wpath1 = VAR.Wstar + (PAR.alpha * VAR.Wstar/(1-PAR.alpha)) * ...
    ( (1 - PAR.alpha * PAR.phi ) * (Kpath1 - VAR.Kstar)/VAR.Kstar ...
      + (PAR.phi - 1) * (Cpath1 - VAR.Cstar)/VAR.Cstar  );
Wpath2 = VAR.Wstar + (PAR.alpha * VAR.Wstar/(1-PAR.alpha)) * ...
    ( (1 - PAR.alpha * PAR.phi ) * (Kpath2 - VAR.Kstar)/VAR.Kstar ...
      + (PAR.phi - 1) * (Cpath2 - VAR.Cstar)/VAR.Cstar  );
Wpath3 = VAR.Wstar + (PAR.alpha * VAR.Wstar/(1-PAR.alpha)) * ...
    ( (1 - PAR.alpha * PAR.phi ) * (Kpath3 - VAR.Kstar)/VAR.Kstar ...
      + (PAR.phi - 1) * (Cpath3 - VAR.Cstar)/VAR.Cstar  );
  
figure(7)
hold on
box off
%ylim([0,1.2])
% xlabel('postshock time (quarters)')
% ylabel('(w(t)-w^*)/dG')
plot(time,(EXO.G/VAR.Wstar)*(Wpath1-VAR.Wstar)/dG,'-','LineWidth',3)
plot(time,(EXO.G/VAR.Wstar)*(Wpath2-VAR.Wstar)/dG,'--','LineWidth',3)
plot(time,(EXO.G/VAR.Wstar)*(Wpath3-VAR.Wstar)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Wpath') 
hold off
pause(1)
  
Rpath1 = PAR.rho + (PAR.rho + PAR.delta) * ...
    ( (Ypath1 - VAR.Ystar)/VAR.Ystar ...
    - (Kpath1 - VAR.Kstar)/VAR.Kstar  );  
Rpath2 = PAR.rho + (PAR.rho + PAR.delta) * ...
    ( (Ypath2 - VAR.Ystar)/VAR.Ystar ...
    - (Kpath2 - VAR.Kstar)/VAR.Kstar  );  
Rpath3 = PAR.rho + (PAR.rho + PAR.delta) * ...
    ( (Ypath3 - VAR.Ystar)/VAR.Ystar ...
    - (Kpath3 - VAR.Kstar)/VAR.Kstar  );

figure(8)
hold on
box off
%ylim([0,1.2])
% xlabel('postshock time (quarters)')
% ylabel('(r(t)-\rho)/dG')
plot(time,(EXO.G/VAR.Rstar)*(Rpath1-PAR.rho)/dG,'-','LineWidth',3)
plot(time,(EXO.G/VAR.Rstar)*(Rpath2-PAR.rho)/dG,'--','LineWidth',3)
plot(time,(EXO.G/VAR.Rstar)*(Rpath3-PAR.rho)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Rpath') 
hold off
pause(1)

figure(9)
hold on
box off
%ylim([0,0.1])
% xlabel('postshock time (quarters)')
% ylabel('\gamma_{C}(t)')
plot(time,100*Cdotpath1./Cpath1,'-','LineWidth',3)
plot(time,100*Cdotpath2./Cpath2,'--','LineWidth',3)
plot(time,100*Cdotpath3./Cpath3,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_gammaCpath') 
hold off
pause(1)

figure(10)
hold on
box off
%ylim([0,0.1])
% xlabel('postshock time (quarters)')
% ylabel('\gamma_{K}(t)')
plot(time,100*Kdotpath1./Kpath1,'-','LineWidth',3)
plot(time,100*Kdotpath2./Kpath2,'--','LineWidth',3)
plot(time,100*Kdotpath3./Kpath3,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_gammaKpath') 
hold off
pause(1)

figure(11)
hold on
box off
%ylim([0,0.1])
% xlabel('postshock time (quarters)')
% ylabel('\gamma_{Y}(t)')
plot(time,100*Ydotpath1./Ypath1,'-','LineWidth',3)
plot(time,100*Ydotpath2./Ypath2,'--','LineWidth',3)
plot(time,100*Ydotpath3./Ypath3,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_gammaYpath') 
hold off
pause(1)

Ipath1 = Ypath1 - Cpath1 - EXO.G - Gpath1;
Ipath2 = Ypath2 - Cpath2 - EXO.G - Gpath2;
Ipath3 = Ypath3 - Cpath3 - EXO.G - Gpath3;

figure(12)
hold on
box off
%ylim([0,1.2])
% xlabel('postshock time (quarters)')
% ylabel('(I(t)-I^*)/dG')
plot(time,(Ipath1-VAR.Istar)/dG,'-','LineWidth',3)
plot(time,(Ipath2-VAR.Istar)/dG,'--','LineWidth',3)
plot(time,(Ipath3-VAR.Istar)/dG,':','LineWidth',3)
legend('\xi_G = 0','\xi_G = 0.02','\xi_G = 0.10','location','best')
% plot([0 1],[0.5 0.5],':','LineWidth',1)
% plot([0.343 0.343],[0 0.5],':','LineWidth',1)
print('-deps2','-loose','F13_08_Ipath') 
hold off
pause(1)

% Impact effects linearized model
dYzero1dG      = (Ypath1(1) - VAR.Ystar)/dG ;
dCzero1dG      = (Czero1 - VAR.Cstar)/dG ;
dIzero1dG      = (Ipath1(1) - VAR.Istar)/dG ;
ElastLGzero1   = ((Lpath1(1) - VAR.Lstar)/VAR.Lstar) / (dG/EXO.G)
ElastRGzero1   = ((Rpath1(1) - VAR.Rstar)/VAR.Rstar) / (dG/EXO.G)
ElastWGzero1   = ((Wpath1(1) - VAR.Wstar)/VAR.Wstar) / (dG/EXO.G)

