clear all
close all
set(0,'defaultlinelinewidth',1.5);

% load data
d = dbload('us_quarterly.csv');

%% get ranges
first_date = get(d.gdp,'first');
last_date  = get(d.gdp,'last');
rng_all    = get(d.gdp,'range');
rng        = qq(1970,1):enddate(d.gdp);  % selected range

%%hezci zobrazeni 
% dat2str(rng)'

%%vyber casove rady
% d.gdp(rng)   % cisla
% d.gdp{rng}   % casova rada

%% take log
% Why?  difference is percentage change

d.l_gdp = 100*log(d.gdp);

% plot graphs
figure
plot(d.gdp)
title('GDP, level, billions of dollars')
grid

figure
plot(d.l_gdp)
title('GDP, 100*log')
grid

figure
plot(rng,d.l_gdp)
title('GDP, 100*log')
grid

%% linerani trend 

d.l_gdp_ltrend = trend(d.l_gdp,'connect=',false);  % iris function
d.l_gdp_lgap   = d.l_gdp - d.l_gdp_ltrend;


figure
plot([d.l_gdp d.l_gdp_ltrend])
legend('log HDP','linearni trend','Location','NorthWest')
title('Linearni trend')
grid

figure
plot(d.l_gdp_lgap)
title('Odchylka od linearniho trendu')
ylabel('%')
grid

%% tempa rustu

% mezictvrtletni, annualizovane
d.dla_gdp  = diff(d.l_gdp)*4;
d.dla_gdp2 = (d.l_gdp - d.l_gdp{-1})*4;   % alternative way

% mezirocni 
d.d4l_gdp  = diff(d.l_gdp,-4);
d.d4l_gdp2 = d.l_gdp - d.l_gdp{-4};       % alternative way


figure
plot(d.dla_gdp)
title('Mezictvrtlentí tempo rustu, (q-on-q), % annualizovane')
grid

figure
plot(d.d4l_gdp)
title('Mezirocni tempo rustu, (y-on-y), %')
grid

% Porovnani s tempem rustu pocitanem z nelog hodnot, fce: pct
d.d4p_gdp  = pct(d.gdp,-4);
d.d4p_gdp2 = (d.gdp / d.gdp{-4} - 1)*100;


figure
plot([d.d4l_gdp d.d4p_gdp])
title('Mezirocni tempo rust, (y-on-y), %')
legend('log approx','not log')
grid

%% Hodrick-Prescott filtr

[d.l_gdp_trend, d.l_gdp_gap] = hpf(d.l_gdp,'lambda=',1600);
% try different lambda

figure
plot([d.l_gdp d.l_gdp_trend])
title('Real GDP, 100*log')
legend('level','trend','Location','NorthWest')
grid

figure
plot([d.l_gdp_gap])
title('Real GDP gap, %')
ylabel('% deviation from trend')
grid

%% Band-pass filter

% frequencies of business cycle  (lower is trend, higher are statistical
% discrepancies)

low  = 8;
high = 32;

[d.l_gdp_bpgap] = bpass(d.l_gdp,[low high]);


figure
plot([d.l_gdp_bpgap d.l_gdp_gap])
title('Real GDP gap, %')
ylabel('% deviation from trend')
legend('Band pass filter','HP filter')
grid


set(gcf,'Papersize',[32 26]);
set(gcf,'PaperPosition',[0 0 32 26]);
saveas(gcf,'filters_comparison','pdf');
%%  statisticke charatekristiky

% tempo rustu
gr.mean = mean(d.d4l_gdp);
gr.std  = std(d.d4l_gdp);
gr.min  = min(d.d4l_gdp);
gr.max  = max(d.d4l_gdp);
gr.ac   = acfn(d.d4l_gdp,5);

gr.pct_pos  = length(find(d.d4l_gdp >0))/length(d.d4l_gdp);   
gr.pct_high = length(find(d.d4l_gdp > gr.mean))/length(d.d4l_gdp);

% hp filtr
hp.mean = mean(d.l_gdp_gap);
hp.std  = std(d.l_gdp_gap);
hp.min  = min(d.l_gdp_gap);
hp.max  = max(d.l_gdp_gap);
hp.ac   = acfn(d.l_gdp_gap,5);

hp.pct_pos  = length(find(d.l_gdp_gap >0))/length(d.l_gdp_gap);   
hp.pct_high = length(find(d.l_gdp_gap > hp.mean))/length(d.l_gdp_gap);

disp('Tempo rustu')
gr
disp('HP filtr')
hp