clear all; 
close all;

sdate = qq(2017,1);
edate = sdate + 15;
rng   = sdate:edate;
frng  = sdate-3:edate;

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

% solve the model
m = solve(m);
m = sstate(m);
get(m,'sstate') 

%% simulate

db = sstatedb(m,sdate-3:sdate);
db.eps_y(sdate) = 1;
% db.eps_pie(sdate) = 1;
%db.eps_i(sdate) = 1;

sdb = simulate(m,db,rng,'deviation', false, 'anticipate',false);
s   = dbextend(db,sdb);

% plot IRFs
f1=figure;  
plot(frng,[s.y_gap s.pie s.i],'linewidth',2,'marker','o');
legend('y','\pi','i');
ylabel('% odchylka od ss');
title('Poptavkovy sok');
grid on;

%% change of inflation target

% create copy of previous model
m2 = m;

% change inflation target
m2.target = 4;

% solve the model and calculte new steady state 
m2 = solve(m2);
m2 = sstate(m2);

% compare steady states of both models
ss  = get(m,'sstatelevel');
ss1 = get(m2,'sstatelevel');
ss&ss1

% simulate the disinflation process starting from original steady state
% set up time range for simulation 
sdate = qq(2017,1);
edate = sdate + 30;

% create steady state database for model with initial target
sdb = sstatedb(m,sdate-3:sdate);

% imulate new model (with new target) but starting from initial database 
sim_db = simulate(m2,sdb,sdate:edate);
s = dbextend(sdb,sim_db);

% plot and interpret the results
dbplot(s,{'pie','y_gap','i','r_gap'},frng,'linewidth',2);

f1=figure;  
plot(frng,[s.y_gap s.pie s.i],'linewidth',2,'marker','o');
legend('y','\pi','i','location','northwest');
ylabel('% odchylka od ss');
title('Zvyseni inflacniho cile');
grid on;

% sacrifice ratio 
sac_ratio = -(cumsum(sim_db.y_gap)/4/(m.target - m2.target));
disp('Sacrifice ratio')
disp([sac_ratio(end)]);