clear all; 
close all;


%% take the model
% 
%  read model file 
m = model('rbc_nonlin_cv8.model');

% calculate analytical steady state
hstar = 0.33;    
kstar = ((1/m.beta - (1-m.delta))/m.alpha)^(1/(m.alpha-1)) * hstar;   
ystar = ((1/m.beta - (1-m.delta))/m.alpha) *kstar;   
cstar = (1- m.alpha*m.delta/(1/m.beta - (1-m.delta))) *ystar;    
rstar = m.alpha*kstar^(m.alpha-1)*hstar^(1-m.alpha); 
wstar = (1-m.alpha)*kstar^(m.alpha)*hstar^(-m.alpha); 
istar = ystar - cstar;

% assign sstate
m.h = hstar;
m.k = kstar;
m.y = ystar;
m.c = cstar;
m.R = rstar;
m.w = wstar;
m.i = istar;
m.z = 0;

% show the sstate
get(m,'sstate')

% let IRIS calculate sstate
m = sstate(m);
get(m,'sstate')

% solve the model
m = solve(m);

%% simulate IRFs

% define interval for simulation
sdate = qq(2000,1);
rng   = sdate:sdate+100;
db    = sstatedb(m,sdate);   % create sstate for all model variables,  
db.epsilon(sdate) = m.sigma;

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

figure;  
plot([s.y s.c s.h s.k s.i],'linewidth',2);
legend('y','c','h','k','i')
grid on;

% relative to ss
figure;  
plot([s.y/m.y s.c/m.c s.h/m.h s.k/m.k s.i/m.i],'linewidth',2);
legend('y','c','h','k','i')
grid on;

%%  simulate shocks in every period

rng2   = sdate:sdate+500;
db2    = sstatedb(m,sdate);   
db2.epsilon = tseries(rng2,@randn)*m.sigma;

s2 = simulate(m,db2,rng,'anticipate',false);

figure;  
plot([s2.y/m.y s2.c/m.c s2.h/m.h s2.k/m.k s2.i/m.i]-1);
legend('y','c','h','k','i')
grid on;