clear all
close all

dm = dbload('us_monthly.csv');
dq = dbload('us_quarterly.csv');

rng_m = mm(1990,1):enddate(dm.unempl);
rng   = qq(1970,1):enddate(dq.gdp);

%% comparison of seasonal adjustment
% unemployment
dm.unempl_sa = x12(dm.unempl_su);

figure
plot(rng_m,[dm.unempl_su dm.unempl_sa])
legend('NSA','SA')
title('Unemployment rate, %')
grid

figure
plot(rng_m,[dm.unempl_sa dm.unempl])
legend('Adjusted by IRIS','Adjusted from source')
title('Unemployment rate, %')
grid

% Task: Do the same for CPI

%% convert to quarterly frequency
d.unr = convert(dm.unempl,'q');
d.ir  = convert(dm.int_rate,'q');
d.cpi = convert(dm.cpi,'q');

% merge two databases
d = dbmerge(d,dq);

%% adjustment 
% get names of the variables
varlist = fieldnames(d);

% take log 
d.l_gdp    = 100*log(d.gdp);
d.l_cons   = 100*log(d.cons);
d.l_invest = 100*log(d.invest);
d.l_real_wage = 100*log(d.wage/d.cpi);
d.l_cpi    = 100*log(d.cpi);
d.pie      = d.l_cpi - d.l_cpi{-4};
d.rir      = d.ir - d.pie{+4};
d.nx_to_gdp  = (d.nx / d.gdp)*100;         % share of NX to GDP, in pct

%% Estimate HP trend and gap for some variables 

[d.l_gdp_trend, d.l_gdp_gap]       = hpf(d.l_gdp,'lambda=',1600);
[d.l_cons_trend, d.l_cons_gap]     = hpf(d.l_cons,'lambda=',1600);
[d.l_invest_trend, d.l_invest_gap] = hpf(d.l_invest,'lambda=',1600);
[d.unr_trend, d.unr_gap]           = hpf(d.unr,'lambda=',1600);

%% calculate correlation coefficient between GDP and other variables (separate function corrfn)

[cc.unr test.unr]        = corrfn(d.l_gdp_gap, d.unr_gap,5);
[cc.cons test.cons]      = corrfn(d.l_gdp_gap, d.l_cons_gap,5);
[cc.invest test.invest]  = corrfn(d.l_gdp_gap, d.l_invest_gap,5);


%% calculate volatilities (STDs)

sd.gdp    = std(d.l_gdp_gap);
sd.unr    = std(d.unr_gap);
sd.cons   = std(d.l_cons_gap);
sd.invest = std(d.l_invest_gap);

% relative to GDP
sd.gdp_to_gdp    = sd.gdp/sd.gdp;
sd.unr_to_gdp    = sd.unr/sd.gdp;
sd.cons_to_gdp   = sd.cons/sd.gdp;
sd.invest_to_gdp = sd.invest/sd.gdp;

%% GDP: trend and gap
figure
plot([d.l_gdp d.l_gdp_trend])
title('Real GDP, 100*log')
legend('level','trend')
grid

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

%%  plot variables and cross-corelation

% unemployment
figure
plot(rng,[d.l_gdp_gap d.unr_gap]);
title('GAPs')
ylabel('% deviation from trend');
legend('GDP','UNR');
grid

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

ps  = cc.unr(:,2);  % phase shift

figure
bar(ps,cc.unr(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('GDP_t vs. UNR_{t+k}')


%% consumption
figure
plot([d.l_gdp_gap d.l_cons_gap]);
title('Consumption')
legend('GDP','C')
ylabel('% deviation from trend');
grid

figure
bar(ps,cc.cons(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('GDP_t vs. C_{t+k}')

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

%% investment

figure
plot([d.l_gdp_gap d.l_invest_gap]);
title('Investment')
legend('GDP','I')
ylabel('% deviation from trend');
grid

figure
bar(ps,cc.invest(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('GDP_t vs. I_{t+k}')

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

%%  standard deviations
disp('Standard deviations')
disp('     GDP        UNR        C        I')
disp([sd.gdp sd.unr  sd.cons sd.invest])
disp('SD relative to GDP')
disp([sd.gdp_to_gdp sd.unr_to_gdp sd.cons_to_gdp sd.invest_to_gdp])