%
%  Program 15.1: Features of an actual demography
%
% Ben J. Heijdra, Foundations of Modern Macroeconomics, Third Ed.
% Oxford University Press
%
% Chapter 15
%
% Ben J. Heijdra
% Groningen, October 26, 2015


clear
clf
clc
close all

%  Load the demographic data
%
%  Source: Human Mortality Database, Cohort 1960, The Netherlands
%  Column 1: Economic age (= biological age - 18 years)
%  Column 2: Survivers (out of 100000 born in 1960)
%  Column 3: Surviving fraction 

DemoData      = importdata('demodata.txt');
EXO.age       = DemoData(:,1);
EXO.survnmbr  = DemoData(:,2);
EXO.survfrcn  = DemoData(:,3);


% Estimate the Boucekkine et al. demography
% Initial guess:

eta0     = 122;
eta1     = 0.068;
x0       = [eta0 eta1];

[PAR.eta_0 PAR.eta_1] = Function15_01_est(EXO,x0);

% Maximu age and predicted surviving fraction

PAR.Dbar         = log(PAR.eta_0) / PAR.eta_1; 
EXO.survfrcnhat  = (EXO.age<PAR.Dbar) .* (PAR.eta_0 - exp(PAR.eta_1 * EXO.age) ) ...
    / (PAR.eta_0 - 1); 
EXO.insmorthat   = (EXO.age<PAR.Dbar) .* (PAR.eta_1 * exp(PAR.eta_1 * EXO.age) ) ...
    ./ (PAR.eta_0 - exp(PAR.eta_1 * EXO.age)); 

EXO.expremlt     = Function15_01_Delta(EXO.age,0,PAR,EXO);  


figure(1)
plot(EXO.age+18,EXO.survfrcn,'-','LineWidth',3)
hold on;
plot(EXO.age+18,EXO.survfrcnhat,':','LineWidth',3)
hold off;
legend('Data','Estimates','Location','Best')
box off
xlim ([18 PAR.Dbar+35])
% ylim ([0 1])
xlabel('biological age (u+18)')
%ylabel('e^{-M(u)}')
print('-deps2','-loose','F15_03')
pause(1)

figure(2)
plot(EXO.age(EXO.age<PAR.Dbar)+18,EXO.insmorthat(EXO.age<PAR.Dbar),'-','LineWidth',3)
hold on;
hold off;
box off
xlim ([18 PAR.Dbar+17])
% ylim ([0 1.2])
xlabel('biological age (u+18)')
%ylabel('\mu(u)')
print('-deps2','-loose','F15_04A')
pause(1)

figure(3)
y_A = EXO.expremlt(30);
y_B = EXO.expremlt(60);
plot(EXO.age+18,EXO.expremlt,'-','LineWidth',3)
hold on;
plot([0 30+17],[y_A y_A],':','LineWidth',3)
plot([30+17 30+17],[0 y_A],':','LineWidth',3)
text(30+19,y_A+2,'A','FontSize',16)
plot([0 60+17],[y_B y_B],':','LineWidth',3)
plot([60+17 60+17],[0 y_B],':','LineWidth',3)
text(60+19,y_B+2,'B','FontSize',16)
hold off;
box off
xlim ([18 PAR.Dbar+17])
% ylim ([0 1.2])
xlabel('biological age (u+18)')
% ylabel('Expected remaining lifetime (\Delta(u,0))')
print('-deps2','-loose','F15_04B')
pause(1)


% Redefine age with fined grid for consumption and asset profiles
EXO.age = linspace(0,100,10001) ;

% Structural parameters

PAR.rho  = 0.04;
EXO.r    = 0.05;
EXO.w    = 1.00;

% Delta terms

EXO.Delta0r      = Function15_01_Delta(0,EXO.r,PAR,EXO);
EXO.Delta0rho    = Function15_01_Delta(0,PAR.rho,PAR,EXO);
EXO.Deltaager    = Function15_01_Delta(EXO.age,EXO.r,PAR,EXO);
EXO.Deltaagerho  = Function15_01_Delta(EXO.age,PAR.rho,PAR,EXO);
VAR.C    = EXO.w * (EXO.Delta0r ./ EXO.Delta0rho) .* exp((EXO.r-PAR.rho) * EXO.age);
VAR.A    = EXO.w * (EXO.Deltaagerho .* VAR.C - EXO.Deltaager) ;

% figure(4)
% plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.C(EXO.age<PAR.Dbar),'-','LineWidth',3)
% hold on;
% hold off;
% box off
% xlim ([18 PAR.Dbar+18])
% % ylim ([0 1.2])
% xlabel('biological age (u+18)')
% %ylabel('\mu(u)')
% print('-deps2','-loose','F15_XX')
% pause(1)
% 
% figure(5)
% plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.A(EXO.age<PAR.Dbar),'-','LineWidth',3)
% hold on;
% hold off;
% box off
% xlim ([18 PAR.Dbar+18])
% % ylim ([0 1.2])
% xlabel('biological age (u+18)')
% %ylabel('\mu(u)')
% print('-deps2','-loose','F15_XX')
% pause(1)

% Solve for u*: the age at which borrowing constraint becomes % binding

ustar0     = 68 ;
fun1       = @(x) PAR.eta_0 * (exp(PAR.rho*x) - 1) + (PAR.rho / (PAR.eta_1 - PAR.rho)) ...
             * (exp(PAR.rho*x) - exp(PAR.eta_1*x)) ...
             - (PAR.rho/EXO.r) * (exp(EXO.r*x) - 1) * (PAR.eta_0 - exp(PAR.eta_1*x));
VAR.ustar      = fsolve(fun1,ustar0);

% Constrained consumption and asset paths

VAR.Cna0   = EXO.w * (PAR.eta_0 - 1) * exp(-(EXO.r - PAR.rho) * VAR.ustar) ...
    ./ (PAR.eta_0 - exp(PAR.eta_1 * VAR.ustar));

VAR.Cna     = (VAR.Cna0 / (PAR.eta_0 - 1)) * (PAR.eta_0 * exp((EXO.r - PAR.rho) * EXO.age) ...
    - exp((EXO.r - PAR.rho + PAR.eta_1) * EXO.age));
VAR.Cna     = (EXO.age<VAR.ustar) .* VAR.Cna + (EXO.age>=VAR.ustar) * EXO.w ;

GAMMA   = (1/(PAR.eta_0 - 1)) * ( (PAR.eta_0/PAR.rho) *(1 - exp(-PAR.rho * EXO.age)) ...
    + (1/(PAR.eta_1 - PAR.rho)) * (1 - exp((PAR.eta_1-PAR.rho) * EXO.age)) );

VAR.Ana     = (EXO.w/EXO.r) * (exp(EXO.r * EXO.age) - 1) ...
    - VAR.Cna0 .* GAMMA .* exp(EXO.r * EXO.age) ;
VAR.Ana     = (EXO.age<VAR.ustar) .* VAR.Ana ;

% figure(6)
% plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Cna(EXO.age<PAR.Dbar),'-','LineWidth',3)
% hold on;
% hold off;
% box off
% xlim ([18 PAR.Dbar+18])
% % ylim ([0 1.2])
% xlabel('biological age (u+18)')
% %ylabel('\mu(u)')
% print('-deps2','-loose','F15_XX')
% pause(1)
% 
% figure(7)
% plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Ana(EXO.age<PAR.Dbar),'-','LineWidth',3)
% hold on;
% hold off;
% box off
% xlim ([18 PAR.Dbar+18])
% % ylim ([0 1.2])
% xlabel('biological age (u+18)')
% %ylabel('\mu(u)')
% print('-deps2','-loose','F15_XX')
% pause(1)

% Constrained consumption and asset paths (ignoring BC)

GAMMAi  = (1/(PAR.eta_0 - 1)) * ( (PAR.eta_0/PAR.rho) *(1 - exp(-PAR.rho * PAR.Dbar)) ...
    + (1/(PAR.eta_1 - PAR.rho)) * (1 - exp((PAR.eta_1-PAR.rho) * PAR.Dbar)) );
VAR.Cnai0  = (1/GAMMAi) * (EXO.w / EXO.r) *(1 - exp(-EXO.r*PAR.Dbar) ) ;
VAR.Cnai   = (VAR.Cnai0 / (PAR.eta_0 - 1)) * (PAR.eta_0 * exp((EXO.r - PAR.rho) * EXO.age) ...
    - exp((EXO.r - PAR.rho + PAR.eta_1) * EXO.age));
VAR.Anai   = (EXO.w/EXO.r) * (exp(EXO.r*EXO.age) -1) - VAR.Cnai0 * GAMMA .* exp(EXO.r*EXO.age) ;

% figure(8)
% plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Cnai(EXO.age<PAR.Dbar),'-','LineWidth',3)
% hold on;
% hold off;
% box off
% xlim ([18 PAR.Dbar+18])
% % ylim ([0 1.2])
% xlabel('biological age (u+18)')
% %ylabel('\mu(u)')
% print('-deps2','-loose','F15_XX')
% pause(1)
% 
% figure(9)
% plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Anai(EXO.age<PAR.Dbar),'-','LineWidth',3)
% hold on;
% hold off;
% box off
% xlim ([18 PAR.Dbar+18])
% % ylim ([0 1.2])
% xlabel('biological age (u+18)')
% %ylabel('\mu(u)')
% print('-deps2','-loose','F15_XX')
% pause(1)

% All figures combined for in the book

figure(10)
plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.C(EXO.age<PAR.Dbar),'-','LineWidth',3)
hold on;
plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Cna(EXO.age<PAR.Dbar),'--','LineWidth',3)
plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Cnai(EXO.age<PAR.Dbar),'-.','LineWidth',3)
plot([VAR.ustar+18 VAR.ustar+18],[0 1],':','LineWidth',2)
text(VAR.ustar+18,1.15,'A','FontSize',16)
hold off;
legend('with annuities', 'without annuities','wrong','location','best')
box off
xlim ([18 PAR.Dbar+18])
% ylim ([0 1.2])
xlabel('biological age (u+18)')
%ylabel('\mu(u)')
print('-deps2','-loose','F15_05A')
pause(1)

figure(11)
plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.A(EXO.age<PAR.Dbar),'-','LineWidth',3)
hold on;
plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Ana(EXO.age<PAR.Dbar),'--','LineWidth',3)
plot(EXO.age(EXO.age<PAR.Dbar)+18,VAR.Anai(EXO.age<PAR.Dbar),'-.','LineWidth',3)
plot([VAR.ustar+18 VAR.ustar+18],[-4 0],':','LineWidth',2)
text(VAR.ustar+18,0.5,'A','FontSize',16)
hold off;
legend('with annuities', 'without annuities','wrong','location','best')
box off
xlim ([18 PAR.Dbar+18])
% ylim ([0 1.2])
xlabel('biological age (u+18)')
%ylabel('\mu(u)')
print('-deps2','-loose','F15_05B')
pause(1)