%%%%% Scrap of Matlab code for optimal growth model

clear all;
close all;

%%%%% parameters

alpha = 0.30;    %% capital's share
beta  = 1/1.05;  %% discount factor
delta = 0.05;    %% depreciation rate
sigma = 1;       %% coefficient relative risk aversion


%%%%% steady state

rho   = 1/beta - 1;                %% discount rate

r     = rho+delta;                 %% steady-state mpk

kbar  = (alpha/r)^(1/(1-alpha));   %% steady-state k
ybar  = kbar^alpha             ;   %% steady-state y
cbar  = ybar - delta*kbar;         %% steady-state c

d2f   = alpha*(alpha-1)*kbar^(alpha-2); %% f''(kbar)


fprintf('\n');
fprintf('Steady state values \n');
fprintf('\n');
fprintf('capital/output       = %7.2f \n', kbar/ybar);
fprintf('consumption/output   = %7.2f \n', cbar/ybar);
fprintf('consumption/capital  = %7.2f \n', cbar/kbar);
fprintf('\n');


%%%%% coefficient matrix

AA = [ 1-beta*d2f*cbar/sigma , d2f*kbar/sigma;
        -cbar/kbar           , 1/beta           ];


fprintf('\n');
fprintf('Coefficient matrix  AA = [a_11 , a_12; a_21 , a_22] \n');
fprintf('\n');
fprintf('a_11   = %7.2f \n', AA(1,1));
fprintf('a_12   = %7.2f \n', AA(1,2));
fprintf('a_21   = %7.2f \n', AA(2,1));
fprintf('a_22   = %7.2f \n', AA(2,2));
fprintf('\n');

%%%%% eigenvalues are roots of quadratic in lambda

lambdas = roots([1,-trace(AA),det(AA)]);

%%%%%% choose stable root

lambda  = min(abs(lambdas));

if abs(lambda)>1,
    display('check roots')
end

%%%%%% solution

psi_kk  = lambda;
psi_ck  = (1/beta - psi_kk)*kbar/cbar;

fprintf('\n');
fprintf('Psi = [psi_kk;psi_ck] \n');
fprintf('\n');
fprintf('psi_kk = %7.2f \n', psi_kk);
fprintf('psi_ck = %7.2f \n', psi_ck);
fprintf('\n');



%%%%% transitional dynamics

k0      = -0.2; %% initial k, log-dev from steady state

T       = 51;   %% horizon

kt      = zeros(T,1);
ct      = zeros(T,1);

kt(1)   = k0;
ct(1)   = psi_ck*k0;

for t=1:T-1,
    
    kt(t+1) = psi_kk*kt(t);
    ct(t+1) = psi_ck*kt(t+1);
    
end


time = (0:1:T-1); %% time index

figure(1)
plot(time,kt,'rd-',time,ct,'bd-')
ylabel('log-deviation from steady state')
xlabel('time')
axis([min(time) max(time) k0 0])
legend('capital','consumption','location','SouthEast')


















