clear all
close all

x = load('us.txt');
t = x(:,1);
gdp = x(:,2);
c = x(:,3);
cpi = x(:,4);
ir = x(:,5);
w = x(:,6);
inv = x(:,7);
n = length(t);

lgdp = log(gdp);
lc = log(c);

lcpi = log(cpi);
lw = log(w);
linv = log(inv);

[lgdp_trend,lgdp_gap] = hp(lgdp,1600);
sd.y = std(lgdp_gap);

%% spotreba
[lc_trend,lc_gap] = hp(lc,1600);

figure
plot(t,[lgdp_gap, lc_gap])
legend('gap gdp','gap c')
grid

cc = corrfn(lgdp_gap,lc_gap,5);  % vypocet kroskorelacniho koeficientu

k = (5:-1:-5)';       % vytvoreni vektoru pro fazovy posun
[cc k]

figure
bar(k,cc)
xlabel('k')
ylabel('korelacni koeficient')
title('HDP_t vs. C_{t+k}')

test.c = [cc k corrtest(cc,n)];    % testovani statisticke vyznamnosti koeficientu korelace
sd.c = std(lc_gap);

disp(' korelacni    fazovy    t-stat    t-crit   vyznamnost')
disp(' koeficient   posun                        a/n=1/0 ')
disp(test.c)

disp('smerodatna odchylka')
disp('      y          c')
disp([sd.y   sd.c]*100)
disp('relativne vuci HDP (std c/std y)') 
disp(sd.c/sd.y)


return

%% investice

[linv_trend,linv_gap] = hp(linv,1600);

figure
plot(t,[lgdp_gap, linv_gap])
legend('gap gdp','gap i')

ci = corrfn(lgdp_gap,linv_gap,5);

figure
bar(k,ci)
xlabel('k')
ylabel('korelacni koeficient')
title('HDP_t vs. I_{t+k}')

test.i = [ci k corrtest(ci,n)];
sd.i = std(linv_gap);