function [X]=get_choleskIRF(M,p,S,coefs_var,nIRF)
% INPUT
% M = # of variables in VAR
% p = # of lags in VAR
% coefs_var = matrix of MxM dimensional coefficients [A1 A2 ...  Ap]: NB
% without the intercept term!
% nIRF = number of steps ahead in IRF
% OUTPUT 
% X is a MxMxnIRF matrix of IRF identified with a choleski ordering

    coef_reshaped=coefs_var;
    % now put the coefficients in pxp matrices that are then used to get the MA
    % representation
    coef3D = zeros(M,M,p+1);
    coef3D(:,:,1) = eye(M);
    for k=1:p
        coef3D(:,:,k+1) = -coef_reshaped(:,(k-1)*M+1:k*M);
    end

    w = invertepolynomialmatrix(coef3D,nIRF);
    X = NaN(M,M,nIRF);
    for i=1:nIRF
        X(:,:,i)=w(:,:,i)*S;
    end