%% 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 + 50;
rng   = sdate:edate;

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

% solve the model
m = solve(m);
% get the steady state
get(m,'sstate') 
% let iris calculate its steady state
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);

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

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

% plot the impulse response funciton for two model variables (l_y_gap and pie)

figure;
subplot(211);
plot(sdb.l_y_gap);
ylabel('% deviation from ss');
title('GDP gap');
subplot(212);
plot(sdb.pie);
title('Inflation');
ylabel('% deviation from ss');

%% simulation with different parameter value
% change value of AR parameter Rho1 (persistence of inflation)
m.Rho1 = 0.85;

% solve the model again
m = solve(m);

% simulate the model to the output gap shock
sdb2 = simulate(m,db,rng,'anticipate',false);

% plot the impulse responses and compare to previous exercise
figure;
subplot(211);
plot(sdb2.l_y_gap);
title('GDP gap');
ylabel('% deviation from ss')
subplot(212);
plot(sdb2.pie);
title('Inflation');
ylabel('% deviation from ss')

%% 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
db2.shock_l_y_gap = tseries(rng,@randn);

% simulate the model
sdb3 = simulate(m,db2,rng,'anticipate',false);

% plot the results (simulated time series)
figure;
subplot(211);
plot(sdb3.l_y_gap);
title('GDP gap');
ylabel('% deviation from ss')
subplot(212);
plot(sdb3.pie);
title('Inflation');
ylabel('% deviation from ss')

