clear all;
close all;


% Definovani strukturalnich parametru modelu
% tj. parametry preferenci (domacnosti) a technologii (firmy)
% strukturalni = policy invariant = nezavisle na hopo

% alpha : kapitalovy podil na vystupu
% beta  : diskontni faktor
% delta : mira depreciace
% sigma : parametr rizikove averze, take (inverzni) mezicasova elast.  substituce 
% gamma : uroven technologie

alpha = .35;
beta  = .98;
delta = .025;
sigma = 2;
gamma = 5;


%  Najdete steady state uroven kapitalu jako funkci strukturalnic  parametru
kstar = ((1/beta - 1 + delta)/(alpha*gamma))^(1/(alpha-1));

%  Definovani poctu diskretnich hodnot k 

gk = 101;
k  = linspace(0.95*kstar,1.05*kstar,gk);


%  vypocitejte matici spotreby (gk x gk) 
%  pro vsechny mozne hodnoty  k_t, k_t+1 

for h = 1 : gk
   for i = 1 : gk
        c(h,i) = gamma*k(h)^alpha + (1-delta)*k(h) - k(i);
        if c(h,i) < 0
            c(h,i) = 0;
        end
       % h je citac pro stavovou promenou  k_t
       % i je citac pro ridici promenou k_t+1
   end
end

%  vypocitejte matici uzitku 

if sigma == 1
            u = log(c);
         else
            u = (c.^(1-sigma) - 1)/(1-sigma);            
end


% definujte pocatecni hodnotovou funkci v (vektro nul)
v = zeros(gk,1);

% nastaveni parametru pro cyklus
convcrit = 1E-9;  % konvergencni kriterium (male cislo)
diff = 1;         % arbitrarne zvolena cislo (vetsi jak convcrit)
iter = 0;         % citac iteraci

one = ones(gk,1); 

while diff > convcrit
    %   pro kazde k_t najdi hodnotu k_t+1, ktera maximalizuje sumu okamziteho
    %   uzitku a diskonovaneho budouciho uzitku
    Tv =  max((u + beta*one*v'),[],2);
    iter = iter + 1;
    diff = abs(Tv - v);
    v = Tv;
end

disp(iter)

% najdi, ktery prvek vektoru  k_t+1 dava optimalni hodnotu 
[Tv,gridpoint] =  max((u + beta*one*v'),[],2);

% najdi rozhodovaci pravidlo pro k_t+1 jako funkce stavu k 
% (prirad pvek vektoru hodnote k_t+1)
k_dr = k(gridpoint);      

figure
plot(k,k_dr,k,k);
xlabel('k_t')
ylabel('k_{t+1}')
title('rozhodovaci pravidlo pro k')


% vypocitejte optimalni rozhodovaci pravidlo jako funkce stavu k_t

c_dr = gamma*k.^alpha + (1-delta)*k - k_dr;


figure
plot(k,c_dr)
xlabel('k_t')
ylabel('c_t')
title('rozhodovaci pravidlo pro c')