% 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);