%
%  Program 16.1: Tragedy of annuitisation
%
% Ben J. Heijdra, Foundations of Modern Macroeconomics, Third Ed.
% Oxford University Press
%
% Chapter 16
%
% Ben J. Heijdra
% Groningen, October 30, 2015


clear
clf
clc
close all

% Parameters
PAR.pie      = 0.30;                       % Death probability
PAR.LEB      = 2 - PAR.pie;                % Life expectancy at birth
n_a          = 0.01;                       % Annual population growth rate
PAR.n        = (1 + n_a)^40 - 1;
PAR.beta     = ((1 + PAR.n)^2 / (2 + PAR.n - PAR.pie)) - 1;
delta_a      = 0.06;                       % Annual depreciation rate
PAR.delta    = 1 - (1 - delta_a)^40;
PAR.alpha    = 0.3;                        % Efficiency parameter capital

% Targets
rstar_a0     = 0.04;                       % Annual interest rate
rstar0       = (1 + rstar_a0)^40 - 1;
ystar0       = 1;                          % Normalized output
kstar0       = PAR.alpha / (rstar0 + PAR.delta);
wstar0       = 1 - PAR.alpha;

% Calibration parameters
PAR.Z_0      = kstar0^(-PAR.alpha);
PAR.rho      = (wstar0 + PAR.pie*(1+rstar0)*kstar0)*(1-PAR.pie) / ((1+PAR.n)*kstar0) ...
    + PAR.pie - 2 ;
rho_a        = (1+PAR.rho)^(1/40)-1;
PAR.s        = (1-PAR.pie)/(2+PAR.rho-PAR.pie); % Savings rate 

% Check that correct steady state is found
StStTY  = @(x) (1+PAR.n) * x - PAR.s * ( (1-PAR.alpha)*PAR.Z_0*x^PAR.alpha ...
    + PAR.pie * (1 - PAR.delta + PAR.alpha * PAR.Z_0 * x^(PAR.alpha-1))* x);
VAR.kbase    = fsolve(StStTY,0.05,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-6,'TolFun',1.0e-6));

VAR.ybase         = PAR.Z_0 * VAR.kbase^PAR.alpha ;
VAR.wbase         = (1-PAR.alpha)*VAR.ybase;
VAR.rbase         = PAR.alpha * VAR.ybase / VAR.kbase - PAR.delta;
VAR.ZYbase        = PAR.pie * (1 + VAR.rbase) * VAR.kbase ;
VAR.CYbase        = (1-PAR.s) *(VAR.wbase + VAR.ZYbase);
VAR.CObase        = (1+VAR.rbase) * PAR.s *(VAR.wbase + VAR.ZYbase);
VAR.LambdaYbase   = log(VAR.CYbase) + ((1 - PAR.pie)/(1 + PAR.rho)) * log(VAR.CObase);

% Model parameters
    disp(' ')
    disp('Model parameters')
    disp(sprintf('      birth rate:                                    %8.4f',PAR.beta))
    disp(sprintf('      death rate:                                    %8.4f',PAR.pie))
    disp(sprintf('      annual population growth rate (percent):       %8.4f',100*n_a))
    disp(sprintf('      annual time preference (percent):              %8.4f',100*rho_a))
    disp(sprintf('      capital share parameter:                       %8.4f',PAR.alpha))
    disp(sprintf('      production function constant:                  %8.4f',PAR.Z_0))
    disp(sprintf('      annual capital depreciation rate (percent):    %8.4f',100*delta_a))
    
% Implied calibration
    disp(' ')
    disp('Implied model parameters and steady state')
    disp(sprintf('      interest factor (1+r):                         %8.4f',VAR.rbase))
    disp(sprintf('      annual interest rate (percent):                %8.4f',100*rstar_a0))
    disp(sprintf('      capital intensity:                             %8.4f',VAR.kbase))
    disp(sprintf('      output:                                        %8.4f',VAR.ybase))
    disp(sprintf('      wage rate:                                     %8.4f',VAR.wbase))
    disp(sprintf('      k/y ratio at annual rate:                      %8.4f',PAR.alpha/(rstar_a0+delta_a)))

% Introduction of perfect annuities

Nperiods       = 20;
time           = linspace(0,Nperiods,Nperiods+1);

% Initialize values under annuitization
VAR.kannu       = zeros(numel(time),1);
VAR.LambdaYannu = zeros(numel(time),1);
VAR.CYannu      = zeros(numel(time),1);
VAR.COannu      = zeros(numel(time),1);
VAR.rannu       = zeros(numel(time),1);
VAR.wannu       = zeros(numel(time),1);
VAR.ZYannu      = zeros(numel(time),1);

% Initial conditions
VAR.kannu(1)    = VAR.kbase;
VAR.wannu(1)    = VAR.wbase;
VAR.ZYannu(1)   = VAR.ZYbase;
VAR.COannu(1)   = VAR.CObase;


% Iterate the fundamental difference equationunder annuitization

for j=2:numel(time)
    
VAR.kannu(j)         = PAR.s/(1+PAR.n) * ( (1-PAR.alpha)*PAR.Z_0 * VAR.kannu(j-1)^PAR.alpha ...
                        + VAR.ZYannu(j-1) );
VAR.wannu(j)         = (1-PAR.alpha) * PAR.Z_0 * VAR.kannu(j)^(PAR.alpha);
VAR.rannu(j)         = PAR.alpha * PAR.Z_0 * VAR.kannu(j)^(PAR.alpha-1) - PAR.delta;
VAR.CYannu(j-1)      = (1-PAR.s) * (VAR.wannu(j-1) + VAR.ZYannu(j-1));
VAR.COannu(j)        = ((1+VAR.rannu(j)) / (1-PAR.pie)) * PAR.s * (VAR.wannu(j-1) + VAR.ZYannu(j-1)) ;
VAR.LambdaYannu(j-1) = log(VAR.CYannu(j-1)) + ((1 - PAR.pie)/(1 + PAR.rho)) ...
                       * log(VAR.COannu(j));
end

figure(1)
plot(time,VAR.kannu,'-b.','MarkerSize',20,'LineWidth',1.5) 
xlim ([0 5])
set(gca,'XTick',0:1:5)
box off
hold on; 
plot(time,VAR.kbase + 0 * time,'--b.','MarkerSize',20,'LineWidth',1.5) 
% plot(time,CapitalPerfectRP,'--b.','MarkerSize',20,'LineWidth',3)  
hold off;
legend('TY','PA','Location','Best')
xlabel('post-shock time  (\tau)')
ylabel('capital intensity (k_{t+\tau})')
print('-deps2','-loose','F16_08_K')
pause(1)

figure(2)
plot(time,VAR.LambdaYannu,'-b.','MarkerSize',20,'LineWidth',1.5) 
xlim ([0 5])
set(gca,'XTick',0:1:5)
box off
hold on; 
plot(time,VAR.LambdaYbase + 0 * time,'--b.','MarkerSize',20,'LineWidth',1.5) 
% plot(time,CapitalPerfectRP,'--b.','MarkerSize',20,'LineWidth',3)  
hold off;
legend('TY','PA','Location','Best')
xlabel('post-shock time  (\tau)')
ylabel('Expected lifetime utility (\Lambda^{Y}_{t+\tau})')
print('-deps2','-loose','F16_08_Lambda')
pause(1)


figure(3)
plot(time,VAR.CYannu,'-b.','MarkerSize',20,'LineWidth',1.5) 
xlim ([0 5])
set(gca,'XTick',0:1:5)
box off
hold on; 
plot(time,VAR.CYbase + 0 * time,'--b.','MarkerSize',20,'LineWidth',1.5) 
% plot(time,CapitalPerfectRP,'--b.','MarkerSize',20,'LineWidth',3)  
hold off;
legend('TY','PA','Location','Best')
xlabel('post-shock time  (\tau)')
ylabel('Consumption during youth (C^{Y}_{t+\tau})')
print('-deps2','-loose','F16_08_CY')
pause(1)

figure(4)
plot(time,VAR.COannu,'-b.','MarkerSize',20,'LineWidth',1.5) 
xlim ([0 5])
set(gca,'XTick',0:1:5)
box off
hold on; 
plot(time,VAR.CObase + 0 * time,'--b.','MarkerSize',20,'LineWidth',1.5) 
legend('TY','PA','Location','Best')
hold off;
% legend('TY','PA','PARP','Location','Best')
xlabel('post-shock time  (\tau)')
ylabel('Consumption during old-age (C^{O}_{t+\tau})')
print('-deps2','-loose','F16_08_CO')
pause(1)



