% inversion of a matrix of polynomials in the lag operator
function inverse = invertepolynomialmatrix(poly,nIRF)
% input 
% poly = matrix of order Mx(M*p+1)
% e.g. in a VAR with M variables of order p poly = [A0 A1 A2 ... Ap]
% the first is the impact matrix
% nIRF = desired horizon for IRFs

M = size(poly,1);
p = size(poly,3) - 1;
polyzero = poly(:,:,1);
invpolyzero =inv(polyzero);
for s = 1:p+1,
   newpoly(:,:,s) = invpolyzero*poly(:,:,s);
end;
polynomialmatrix = - newpoly(:,:,2:p+1);
% generate companion matrix A
A = zeros(M*p,M*p);
A(M+1:M*p,1:M*(p-1)) = eye(M*(p-1));
% fill the first M rows of the companion matrix with A0^(-1)A1 A0^(-1)A2 ... A0^(-1)Ap
for j = 1:p
   A(1:M,(j-1)*M+1:j*M) = polynomialmatrix(:,:,j);
end
inverse = zeros(M,M,nIRF);
D = eye(M*p);
for j = 1:nIRF 
   inverse(:,:,j) = D(1:M,1:M)*invpolyzero;
   D = A*D;
end
