%% 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 output gap shock (normally distributed) 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;