
clear all
close all

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

lgdp = log(gdp);
lc = log(c);
lcpi = log(cpi);
lw = log(w);
linv = log(inv);

[lgdp_trend,lgdp_gap] = hp(lgdp,1600);  % estimation of trend and gap
[lc_trend,lc_gap] = hp(lc,1600);
[linv_trend,linv_gap] = hp(linv,1600);

figure
plot(time,[lgdp_gap, lc_gap]*100)
legend('gdp gap','c gap')
ylabel('in %')

figure
plot(time,[lgdp_gap, linv_gap]*100)
legend('gdp gap','inv gap')
ylabel('in %')

c1 = corrfn(lgdp_gap,lc_gap,5);    % cross-correlation coefficient
c2 = corrfn(lgdp_gap,linv_gap,5);

k = (5:-1:-5)';
[k, c1]

figure
bar(k,c1)
title('gdp gap_{t} x c gap_{t+k}')
xlabel('k')
ylabel('corr. coef.')

figure
bar(k,c2)
title('gdp gap_{t} x inv gap_{t+k}')
xlabel('k')
ylabel('corr. coef.')

std_y = std(lgdp_gap);    % standard deviation (mesures volatility of variable)
std_c = std(lc_gap);
std_i = std(linv_gap);

share_cy = c./gdp;       % share on gdp
share_iy = inv./gdp;
avg_cy = mean(share_cy);
avg_iy = mean(share_iy);

disp('Share on GDP:   C   I')
disp([avg_cy  avg_iy]*100)
disp('Standard deviation: GDP  C   I')
disp([std_y  std_c  std_i]*100)

lrw = lw - lcpi;
% rw = w ./cpi; lrw = log(rw);     %alternative way of calculation

[lrw_trend,lrw_gap] = hp(lrw,1600);

c3 = corrfn(lgdp_gap,lrw_gap,5);

figure
plot(time,[lgdp_gap, lrw_gap]*100)
legend('gdp gap','rw gap')
ylabel('in %')

figure
bar(k,c3)
title('gdp gap_{t} x rw gap_{t+k}')
xlabel('k')
ylabel('corr. coef.')

% inflation rate
pie4 = lcpi(5:end) - lcpi(1:end-4);           % y/y inflation rate
pie4 = [nan(4,1);pie4];

[pie4_trend,pie4_gap] = hp(pie4,1600);
c4 = corrfn(lgdp_gap,pie4_gap,5);

figure
plot(time,[lgdp_gap, pie4_gap]*100)
legend('gdp gap','pie gap')
ylabel('in %')

figure
bar(k,c4)
title('gdp gap_{t} x pie gap_{t+k}')
xlabel('k')
ylabel('corr. coef.')

n = length(time);
stats = corrtest(c1,n);

test.c = [c1 k corrtest(c1,n)];
test.i = [c2 k corrtest(c2,n)];
test.rw = [c3 k corrtest(c3,n)];
test.pie = [c4 k corrtest(c4,n-4)]; % minus 4 because of nans

% real interest rate
rir = ir(1:end-4) - pie4(5:end);
rir = [rir;nan(4,1)];

[rir_trend, rir_gap] = hp(rir,1600);
c5 = corrfn(lgdp_gap,rir_gap,5);

figure
plot(time,[lgdp_gap, rir_gap]*100)
legend('gdp gap','rir gap')
ylabel('in %')

figure
bar(k,c5)
title('gdp gap_{t} x rir gap_{t+k}')
xlabel('k')
ylabel('corr. coef.')

test.rir = [c5 k corrtest(c5,n-4)]; % minus 4 because of nans








