function [C,BigA,pi_hat,Y,X,Y_initial,Y_tilde_fit,err, Y_tilde, X_tilde] = wold_irf_LP(vardata,n,p,c,hor,st_f)

    % Reduced-Form VAR estimation
    [pi_hat,Y,X,Y_initial,Y_tilde_fit,err, Y_tilde, X_tilde]=VAR_LP(vardata,p,c,st_f);
    
    % Companion matrix representation:
    BigA=[pi_hat(c+1:end,:)'; eye(n*p-n) zeros(n*p-n,n)]; % np x np matrix
    
    % Wold representation IRFs:
    C=zeros(n,n,hor); % reduced-form IRFs
    
    for j=1:hor
        BigC=BigA^(j-1);
        C(:,:,j)=BigC(1:n,1:n); % IRF of the Wold representation
    end

end