clear all
close all

x = load('us.txt');
lgdp = log(x(:,2));
lcpi = log(x(:,4));
ir = x(:,5)/100;

pie = (lcpi - lag(lcpi,1))*4;   % quarter-on-quarter inflation rate

[lgdp_trend, lgdp_gap] = hp(lgdp, 1600);
[pie_trend, pie_gap] = hp(pie, 1600);
[ir_trend, ir_gap] = hp(ir, 1600);

rir_gap = lag(ir_gap,1) - pie_gap;   

% estimation of parameters of individual equations
X1 = [ lag(lgdp_gap,1)  rir_gap ];
a = ols(lgdp_gap,X1);

X2 = [lag(pie_gap,1), lgdp_gap];
b = ols(pie_gap, X2);

X3 = [lag(ir_gap,1), pie_gap, lgdp_gap];
c = ols(ir_gap, X3);

% structural parameters
alpha = a(1,1);
beta = a(2,1);
gamma = b(1,1);
delta = b(2,1);
omega = c(1,1);
kappa = c(2,1);
lambda = c(3,1);

% matrix form of VAR
A = [1 beta 0
    -delta 1 0
    -lambda -kappa 1];

B = [alpha 0  beta
    0 gamma 0
    0  0  omega];
 

C = inv(A)*B;
D = inv(A);

eigenvalues = abs(eig(C))    % stability of system 

z = zeros([3,20]);  % trajectory of variables (output, inflation, interest rate), 20 periods
epsilon = [1;0;0];  % shock in first period: demand, supply, monetary
 
% simulation of impuslse responses

z(:,1) = D*epsilon;
 
for i = 2:20
     z(:,i) = C * z(:,i-1);
end
 
figure
plot(z')
title('impulse respones')
xlabel('time')
ylabel('deviation from equilibrium')
legend('y gap','\pi gap','ir gap')
grid
 
 
 










