%
%  Program 17.2: Infinitely-lived  agents and the consumption and savings
%
% Ben J. Heijdra, Foundations of Modern Macroeconomics, Third Ed.
% Oxford University Press
%
% Chapter 17
%
% Ben J. Heijdra
% Groningen, November 15, 2015


clear
clf
clc
close all

% Parameters
delta_a      = 0.06;                       % Annual depreciation rate
delta        = 1 - (1 - delta_a)^25;
alpha        = 0.3;                        % Efficiency parameter capital

% Targets
r_a0         = 0.04;                       % Annual interest rate
r            = (1 + r_a0)^25 - 1;
rho_a        = 0.03;
rho          = (1 + rho_a)^25 - 1;
beta         = 1/(1+rho);
ystar0       = 1;                          % Normalized output
kstar0       = alpha / (r + delta);
w            = 1 - alpha;

% Optimal plans -- deterministic

H_1 = w * (1 + 1/(1+r) + 1/(1+r)^2 );
C_1 = H_1 / (1 + beta + beta^2);
C_2 = beta *(1+r)* C_1;
C_3 = beta *(1+r)* C_2;

A_2 = w - C_1;
A_3 = (1+r)*A_2 + w - C_2;
A_4 = (1+r)*A_3 + w - C_3;

% Optimal plans -- stochastic

e_1 = 0.75;
e_2 = 1.00;
e_3 = 1.25;

p_11 = 0.70;
p_12 = 0.25;
p_13 = 1 - p_11 - p_12;
p_21 = 0.10 ;
p_22 = 0.80 ;
p_23 = 1 - p_21 - p_22;
p_33 = 0.70;
p_32 = 0.25;
p_31 = 1 - p_32 - p_33;

pi_11 = 0.00;
pi_12 = 1.00;
pi_13 = 0.00;

pi_21 = p_21;
pi_22 = p_22;
pi_23 = p_23;

pi_31 = p_21 * (p_11 + p_22) + p_23 * p_31;
pi_32 = p_21 * p_12 + p_22^2 + p_23 * p_32 ;
pi_33 = p_21 * p_13 + p_23 * (p_22 + p_33) ;

ExpH_1 = w + (w/(1+r)) * (pi_21 * e_1 + pi_22 * e_2 + pi_23 * e_3) ...
    + w * (1/(1+r))^2 * (pi_31 * e_1 + pi_32 * e_2 + pi_33 * e_3) ;


% Long-run distribution

P = zeros(3,3);
P(1,1) = p_11;
P(1,2) = p_12;
P(1,3) = p_13;
P(2,1) = p_21;
P(2,2) = p_22;
P(2,3) = p_23;
P(3,1) = p_31;
P(3,2) = p_32;
P(3,3) = p_33;

pi1  = [0.00 1.00 0.00]';
pi2  = (pi1' * P)'
pi3  = (pi2' * P)'
pi4  = (pi3' * P)'
pi5  = (pi4' * P)'
pi6  = (pi5' * P)'
pi7  = (pi6' * P)'
pi8  = (pi7' * P)'
pi9  = (pi8' * P)'
pi10 = (pi9' * P)'
pi11 = (pi10' * P)'
pi12 = (pi11' * P)'
pi13 = (pi12' * P)'
pi14 = (pi13' * P)'
pi15 = (pi14' * P)'
pi16 = (pi15' * P)'
pi17 = (pi16' * P)'



% Compute long-run unconditional probabilities

% by hand
Delta = zeros(2,2);
Delta(1,1) = p_11 - 1;
Delta(1,2) = p_21;
Delta(2,1) = p_33 - p_11;
Delta(2,2) = p_22 + p_33 - 2;

pi_bar = Delta^(-1) * [0 ; p_33-1];
pi_1bar = pi_bar(1,1)
pi_2bar = pi_bar(2,1)
pi_3bar = 1 - pi_2bar - pi_1bar

% with Matlab routine 
% V = right eigenvectors in columns
% D = matrix with eigenvalues on main diagonal
% W = left eigenvectors in columns

[V D W] = eig(P);

pie_1bar = W(1,1) / (ones(1,3) * W(:,1))
pie_2bar = W(2,1) / (ones(1,3) * W(:,1))
pie_3bar = W(3,1) / (ones(1,3) * W(:,1))




% Model parameters
%     disp(' ')
%     disp('Exogenous variables')
%     disp(sprintf('      initial assets A_1:                            %8.4f',0))
%     disp(sprintf('      initial human wealth H_1:                      %8.4f',H_1))
%     disp(sprintf('      initial expected human wealth E_1H_1:          %8.4f',ExpH_1))
%     disp(sprintf('      wage w_t:                                      %8.4f',w))
%     disp(sprintf('      interest rate r_t:                             %8.4f',r))
%         
%     disp('Deterministic choices')
%     disp(sprintf('      consumption C_1:                               %8.4f',C_1))
%     disp(sprintf('      consumption C_2:                               %8.4f',C_2))
%     disp(sprintf('      consumption C_3:                               %8.4f',C_3))
%     disp(sprintf('      assets A_2:                                    %8.4f',A_2))
%     disp(sprintf('      assets A_3:                                    %8.4f',A_3))
%     disp(sprintf('      assets A_4:                                    %8.4f',A_4))
%     
%     disp('Stochastic choices - Traditional method')
%     disp(sprintf('      consumption C_1:                               %8.4f',C_1stoch))
%     disp(sprintf('      consumption C_2(e_1):                          %8.4f',C_2stoch_e_1))
%     disp(sprintf('      consumption C_2(e_2):                          %8.4f',C_2stoch_e_2))
%     disp(sprintf('      consumption C_2(e_3):                          %8.4f',C_2stoch_e_3))
%     disp(sprintf('      assets A_2:                                    %8.4f',A_2stoch))
%     disp(sprintf('      assets A_3(e_1):                               %8.4f',A_3stoch_e_1))
%     disp(sprintf('      assets A_3(e_2):                               %8.4f',A_3stoch_e_2))
%     disp(sprintf('      assets A_3(e_3):                               %8.4f',A_3stoch_e_3))
%     disp(sprintf('      consumption C_3(e_1,e_1):                      %8.4f',C_3stoch_e_1_e_1))
%     disp(sprintf('      consumption C_3(e_1,e_2):                      %8.4f',C_3stoch_e_1_e_2))
%     disp(sprintf('      consumption C_3(e_2,e_1):                      %8.4f',C_3stoch_e_2_e_1))
%     disp(sprintf('      consumption C_3(e_2,e_2):                      %8.4f',C_3stoch_e_2_e_2))
%     disp(sprintf('      consumption C_3(e_2,e_3):                      %8.4f',C_3stoch_e_2_e_3))
%     disp(sprintf('      consumption C_3(e_3,e_2):                      %8.4f',C_3stoch_e_3_e_2))
%     disp(sprintf('      consumption C_3(e_3,e_3):                      %8.4f',C_3stoch_e_3_e_3))
% 
%     disp('Stochastic choices - DP method')
%     disp(sprintf('      consumption C_1:                               %8.4f',C_1_e_2(1)))
%     disp(sprintf('      assets A_2:                                    %8.4f',Aplus_1_e_2(1)))
% 
