% zadani cviceni 6

clear all;   % vymaze promenne z pameti
close all;   % zavre vsechna okna
tic          % stopky (start)

alpha   =  0.35;     % kapitalovy podil
beta    =  0.96;     % diskontni faktor
delta   =  0.08;     % mira depreciace
psi     =  1.64;      % parametr pro disutilit z prace

% steady statove hodnoty promennych pred reformou

hstar = 36/105;    % odpracovane hodiny
kstar = ((1/beta - (1-delta))/alpha)^(1/(alpha-1)) * hstar;    % kapital
ystar = ((1/beta - (1-delta))/alpha) *kstar;    % vystup
istar = alpha*delta / (1/beta - (1-delta))*ystar;    % investice
cstar = ystar - istar;    % spotreba
wstar = (1-alpha)*kstar^(alpha)*hstar^(-alpha);    % mzda
rstar = alpha*kstar^(alpha-1)*hstar^(1-alpha);    % najemni cena kapitalu
lstar = 1-hstar;    % volny cas
incstar = hstar*wstar;  % pracovni duchod
ustar = log(cstar)+psi*log(lstar);    % uzitek

disp('kstar ystar istar cstar wstar rstar lstar incstar ustar')
disp([kstar ystar istar cstar wstar rstar lstar incstar ustar])
%% steady statove hodnoty promennych po zavedeni reformy (povinne snizeni
% pracovni doby

hstar2 = 30/105;    % odpracovane hodiny
kstar2 = ((1/beta - (1-delta))/alpha)^(1/(alpha-1)) * hstar2;    % kapital
ystar2 = ((1/beta - (1-delta))/alpha) *kstar2;    % vystup
istar2 = alpha*delta / (1/beta - (1-delta))*ystar2;    % investice
cstar2 = ystar2 - istar2;    % spotreba
wstar2 = (1-alpha)*kstar2^(alpha)*hstar2^(-alpha);    % mzda
rstar2 = alpha*kstar2^(alpha-1)*hstar2^(1-alpha);    % najemni cena kapitalu
lstar2 = 1-hstar2;    % volny cas
incstar2 = hstar2*wstar2;  % pracovni duchod
ustar2 = log(cstar2)+psi*log(lstar2);    % uzitek

disp('kstar2 ystar2 istar2 cstar2 wstar2 rstar2 lstar2 incstar2 ustar2')
disp([kstar2 ystar2 istar2 cstar2 wstar2 rstar2 lstar2 incstar2 ustar2])
lambda = cstar/cstar2*(lstar/lstar2)^psi;    % kompenzacni konstanta
 

%% Nalezeni rozhodovaciho pravidla
%  (pouzijte iteraci hodnotove funkce k nalezeni rozhodovaciho pravidla po
%  zavedeni reformy)


gk  = 101;      % pocet bodu (grid) pro ridici a stavovou promenou k
k = linspace(.7*kstar,1.3*kstar,gk)';  % oblast kapitalu kolem steady statove hodnoty


% matice spotreby a matice uzitku (preddefinovani)
c = zeros(gk,gk);
u = zeros(gk,gk);

warning off; % vypnuti  varovani pri logaritmovani nuly (jen to otravuje)


%  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) = k(h)^alpha*hstar2^(1-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 
u = log(c) + psi*(1-hstar2);
  

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

% nastaveni parametru pro cyklus
convcrit = 1E-9;  % konvergencni kriterium (male cislo)
diff = 1;         % arbitrarne zvolene 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);      
c_dr = k.^alpha*hstar2^(1-alpha) + (1-delta)*k - k_dr;



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

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


return
%% Simulace 

% Predpokladame, ze se ekonomika pred reformou nachazi ve steady statove
% tj. pocatecni hodnota kapitalu je hodnota puvodniho steady statu

ksim(1,1) = kstar;

% delka gridu (vzdalenost mezi sousednimi hodnotami k) 
% vyuzimeme pri simulaci 
grid_length = (k(gk)-k(1))/(gk-1);
 
% cyklus pro simulaci      
for i = 1 : 50
    n = (ksim(i,1)-k(1,1))/grid_length+1;  % nalezeni pozice soucasne hodnoty kapitalu v gridu 
    n = round(n);                          % zaokrouhleni (pozice je vzdy cele cislo)
    ksim(i+1,1) = k_dr(n);                 % prirazeni hodnoty kapitalu k_t+1 dle rozhodovaci pravidla
    csim(i,1) = c_dr(n);               % vypocet spotreby dle rozhodovaci pravidla
    ysim(i,1) = [];                    % vypocet ostatnich simulovanych hodntot
    isim(i,1) = [];               
    rsim(i,1) = [];
    wsim(i,1) = [];
    incsim(i,1) = []; 
    usim(i,1) = [];
end

% Vykresleni prechodu pro makroekonomicke promenne 
% (zacatek je v puvodnim steady statu)

t = 0:50;    % vektor casu pro vykresleni
k_t=[ksim];  % vektor kapitalu je o jeden prvek delsi
c_t=[cstar;csim];
y_t=[];
r_t=[];
w_t=[];
i_t=[];
inc_t = [];
u_t=[];

% Vykresleni




    
