%% start ... houskeeping
clear all; 
close all;


%% set dates for simulation
% 
% specify starting and final date for simulation (impulse responses)

sdate = qq(2016,1);
edate = sdate + 80;
rng   = sdate:edate;

%% take the model
% 
%  read model file 
m = model('rbc2.model','linear=',true);

% solve the model
m = solve(m);
m = sstate(m);
get(m,'sstate')
%% simulate one-time shock
% create initial database with all the variables at their steady-states
db = sstatedb(m,sdate);
db = zerodb(m,sdate);  % alternativni prikaz

% add unit one-off  shock output gap in the second period of simulation range
db.epsilon(sdate+1) = 1;

% simulate reaction of the variables to one-off surprise shock in output gap
sdb = simulate(m,db,rng,'anticipate',false);

% plot the impulse response funciton for two model variables (l_y_gap and pie)
figure;
plot([sdb.y sdb.c sdb.h sdb.k sdb.z sdb.i sdb.w],'linewidth',2);
ylabel('% deviation from ss');
title('Impulzni odezva');
legend('y','c','h','k','z','i','w');
grid;

%% simulation to series of shocks

% create new database as a copy of the previous one
db2 = db;

% now generete series of normally distributed technology shocks for the
% whole period specified for simulation
rng2 = sdate:sdate+500;
db2.epsilon = tseries(rng2,@randn)*m.sigma;

% simulate the model
sdb2 = simulate(m,db2,rng2,'anticipate',false);

% plot the results (simulated time series)
figure;
plot([sdb2.y sdb2.c sdb2.h sdb2.i],'linewidth',2);
title('Simulace');
ylabel('% deviation from ss');
legend('y','c','h','i');
grid on;

%% calculate statistics

% volatility
sd.y    = std(sdb2.y);
sd.c    = std(sdb2.c);
sd.i    = std(sdb2.i);
sd.h    = std(sdb2.h);

% relative to Y
sd.c2y = sd.c/sd.y;
sd.i2y = sd.i/sd.y;
sd.h2y = sd.h/sd.y;

disp('Standard deviation: y c i h')
disp([sd.y sd.c sd.i sd.h])
disp('Standard deviation relative to y: c i h')
disp([sd.c2y sd.i2y sd.h2y])

% correlation
[cc.c test.c]  = corrfn(sdb2.y, sdb2.c,5);
[cc.i test.i]  = corrfn(sdb2.y, sdb2.i,5);
[cc.h test.h]  = corrfn(sdb2.y, sdb2.h,5);
[cc.R test.R]  = corrfn(sdb2.y, sdb2.R,5);
[cc.w test.w]  = corrfn(sdb2.y, sdb2.w,5);
[cc.hw test.hw]  = corrfn(sdb2.h, sdb2.w,5);


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

% figure
% bar(ps,[cc.c(:,1) cc.i(:,1) cc.h(:,1) cc.w(:,1) cc.R(:,1)])
% legend('c','i','h','w','R')
% xlabel('lead                    lag')
% ylabel('corr. coef.')

subplot(221)
bar(ps,cc.i(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('Y vs. I_{t+k}')
subplot(222)
bar(ps,cc.c(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('Y vs. C_{t+k}')
subplot(223)
bar(ps,cc.h(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('Y vs. H_{t+k}')
subplot(224)
bar(ps,cc.R(:,1))
xlabel('lead                    lag')
ylabel('corr. coef.')
title('Y vs. R_{t+k}')


