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
% theta : parametr rizikove averze, take (inverzni) mezicasova elast.  substituce 
% gamma : uroven technologie

alpha = .35;
beta  = .98;
delta = .025;
theta = 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 
g = 3;

%  Mozne stavy kapitalu
k = [2.85
     3.00
     3.15]

% (a)
% Vypocitejte matici spotreby c (3 × 3)   pro vsechny mozne kombinace hodnot k_t a  k_t+1

% vypocet v cyklu
% for i = 1 : g
%    for j = 1 : g
%       c(i,j) = gamma*k(i)^alpha + (1-delta)*k(i) - k(j);
%       % i je citac pro stavovou promenou  k_t
%       % j je citac pro ridici promenou k_t+1
%    end
% end

% efektivni vypocet
c = (gamma*k.^alpha + (1-delta)*k)*ones(1,3) - ones(3,1)*k';

% prirazeni uzitku 
% vypocet v cyklu

% for i = 1 : g
%    for j = 1 : g
%       u(i,j) = (c(i,j)^(1-theta) - 1)/(1-theta);    
%         i je citac pro stavovou promenou  k_t
%         j je citac pro ridici promenou k_t+1
%    end
% end

% efektivni vypocet
u = (c.^(1-theta)-1) / (1-theta);

% (b)
% hodnotova funkce 
v = [167.6
     168.1
     168.6];
 
% Transformace vektoru v na matici 3 x 3, kde kazdy radek je identicky

% definovani vektoru jednicek
one =   [1
         1
         1];

% prenasobeni vektoru jednicek a transponovaneho vektoru v
tmp1 = one*v';

% alternativne pomoci fce repmat
tmp0 = repmat(v',[3,1]);

% Vykresleni vysledku 
disp(''); % blank line
disp('v transformovany:')
disp(tmp1);

% Vypocet souctu {u(k_t,k_t+1) + beta*v(k_t+1)}
tmp2 = u + beta*one*v';

% Vykresleni vysledku
disp(''); % blank line
disp('Matice ve slozenych zavorkach:')
disp(tmp2);

% (c)
%  Nalezeni nejvetsi hodnoty v kazdem radku
%  i.e. pro kazde k_t najdi hodnotu k_t+1, ktera maximalizuje cilovou funkci

Tv = max((u + beta*one*v'),[],2);

disp(''); % blank line
disp('Vysledna hodnota z maximalizace pro kazde k_t:')
disp(Tv);

% (d)
% Nalezeni optimalniho rozhodovaciho pravidla  (zatim jen jeden krok)

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

disp(''); % blank line
disp('Prvek vektoru k_t+1 pri optimalnim vyberu:')
disp(i);

% (optimalni) rozhodovaci pravidlo pro kapital (prevod poradi vektoru na hodnotu kapitalu k_t+1)
kdr = k(i);

disp(''); % blank line
disp('rozhodovaci pravidlo pro k_t+1:')
disp(kdr);

% (optimalni) rozhodovaci pravidlo pro spotrebu, jako funkce stavu k
cdr = gamma*k.^alpha + (1-delta)*k - kdr;


% (e)
% Vyse uvedene v cyklu 

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

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

% rozhodovaci pravidla
kdr = k(i);
cdr = gamma*k.^alpha + (1-delta)*k - kdr;