clear all
close all
start_time = cputime;

%% Initial Setting
beta=0.96;          % individualni discontni faktor   
delta=0.08;         % mira depreciace
sigma=1;            % koeficient relativni averze k riziku
alfa=0.36;          % podil kapitalu v produkcni funkci
 

k_max = 7;
k_min = 0.1;
nk = 500;

k = linspace(k_min,k_max,nk);

    
% produkcni funkce
y = k.^alfa;

% uzitkova funkce
% rozpoctove omezeni: c = y + (1 - delta)k - k_prime
I = ones(size(k));
C(:,:) = I'*(y + (1-delta)*k) - k'*I;



U_C = log(max(C,1e-12));        % podminka nezapornosti C > 0

% Solving the agent's problem

frame=zeros(1,nk);         
TV=ones(size(frame));           % pocatecni odhad pro VF
V=zeros(size(frame));
step=1;
V_norm=1e-9;

figure
%% Value Function Iteration
tic
while and(norm(TV - V)>V_norm, step<5000)
V=TV;

EV = repmat(V',[1,nk]);     % zde zbytecny krok, tohle vyuzijete pri vetsi dimenzi

    RHS = U_C + beta*EV;
    RHS(C<=0) = -1e30;      % jeste jednou aplikace podminky nezapornosti C > 0

    [TV, TV_index] = max(RHS);    % Bellman equation
    k_pol = k(TV_index);

    subplot(2,1,1)
    plot(k,V)
    title(['Hodnotova funkce - iterace ' num2str(step)])
    subplot(2,1,2)
    plot(k,k_pol,k,k)
    title('Rozhodovaci pravidlo - kapital')

    frame = getframe;

step=step+1;
end
toc
step
VF = TV;


%% Policy Functions
% Policy functions
k_pol = k(TV_index);

c_pol = y(TV_index) + (1 - delta)*k - k_pol;

% Obrazky
figure
plot(k,VF)
title('Hodnotova funkce')

figure
plot(k,k_pol,k,k)
title('Rozhodovaci pravidlo - kapital')

figure
plot(k,c_pol,k,k)
title('Rozhodovaci pravidlo - spotreba')

total_time = cputime - start_time