% Calibrated parameter values

clear;
beta=0.99;
sigma=1;
chi=1.55;
eta=0;
theta=2.064;
omega=0.5;
alpha=3;
delta=1.5;
rhov=0.5;
rhou=0.8;
sigmav=1;
sigmau=0.5;

% Calculate kappa 

kappa=(1-omega)*(1-beta*omega)/(alpha*omega);

% Define state space matrices

A0=zeros(4,4);
A0(1,1)=1;
A0(2,2)=1;
A0(3,3)=1;
A0(3,4)=sigma^-1;
A0(4,4)=beta;

A1=zeros(4,4);
A1(1,1)=rhov;
A1(2,2)=rhou;
A1(3,1)=sigma^-1;
A1(3,3)=1;
A1(3,4)=sigma^-1*delta;
A1(4,2)=-1;
A1(4,3)=-kappa;
A1(4,4)=1;

B0=zeros(4,2);
B0(1,1)=1;
B0(2,2)=1;

% Calculate alternative state space matrices

A=inv(A0)*A1;
B=inv(A0)*B0;

% Jordan decomposition of A

[p,lambda]=eig(A);
pstar=inv(p);

% Sort eigenvalues and eigenvectors in ascending order

val=diag(lambda);
t=sortrows([val p'],1);
lambda=diag(t(:,1));
p=t(:,2:5)';
pstar=inv(p);

% Partition matrices

LAMBDA1=lambda(1:2,1:2);
LAMBDA2=lambda(3:4,3:4);

P11=pstar(1:2,1:2);
P12=pstar(1:2,3:4);
P21=pstar(3:4,1:2);
P22=pstar(3:4,3:4);

R=pstar*B;

% Calculate impulse response functions

h=20;
irf_v=zeros(4,h);
irf_u=zeros(4,h);
irf_v(1:2,1)=inv(P11-P12*inv(P22)*P21)*R(1:2,:)*[1;0];
irf_u(1:2,1)=inv(P11-P12*inv(P22)*P21)*R(1:2,:)*[0;1];

i=1;
for i=1:(h-1);
    irf_v(1:2,i+1)=inv(P11-P12*inv(P22)*P21)*LAMBDA1*(P11-P12*inv(P22)*P21)*irf_v(1:2,i);
    irf_u(1:2,i+1)=inv(P11-P12*inv(P22)*P21)*LAMBDA1*(P11-P12*inv(P22)*P21)*irf_u(1:2,i);
end;

irf_u(3:4,:)=-inv(P22)*P21*irf_u(1:2,:);
irf_v(3:4,:)=-inv(P22)*P21*irf_v(1:2,:);

irf_v=real(irf_v);
irf_u=real(irf_u);