function out = irf_fevd(fac_est_out, est_par, hor, tcodevec)
% Compute IRFs and FEVDs. 
% This is adapated from James H. Stock and Mark W. Watson (2016).

% EXTRACT INPUTS
lam_mat       = fac_est_out.lam_mat;
uar_coef_mat  = fac_est_out.uar_coef_mat;
uar_ser_mat   = fac_est_out.uar_ser_mat;
varout        = fac_est_out.varout;
n_uarlag      = est_par.n_uarlag;
[dnobs, nfac] = size(fac_est_out.fac);
n_series      = size(fac_est_out.lam_mat,1);
n_shocks      = size(varout.coef.G,2);

%--- compute IRFs and FEVDs from factor VAR 
coef_y_fac = varout.coef; 
imp_fac_mat = varirf_ss(coef_y_fac,0,hor);    

vcomp_fac_mat = imp_fac_mat.^2;
vtotal_fac_mat = NaN(nfac,hor);

for is = 1:nfac
    tmp = reshape(vcomp_fac_mat(is,:,:),n_shocks,hor);
    tmp = cumsum(tmp,2);
    vcomp_fac_mat(is,:,:) = tmp;
    vtotal_fac_mat(is,:) = sum(tmp);
end

vfrac_fac_comp_mat = NaN(nfac,n_shocks,hor);
for i = 1:n_shocks
    tmp = reshape(vcomp_fac_mat(:,i,:),nfac,hor);
    tmp = tmp./vtotal_fac_mat;
    vfrac_fac_comp_mat(:,i,:) = tmp;
end

%-- compute IRFs for variables to factor shocks
coef_y_fac = varout.coef; 
coef_y_fac.Q = lam_mat*coef_y_fac.Q;
imp_y_fac_mat = varirf_ss(coef_y_fac,0,hor);        % irf to factor shocks

%-- compute IRFs wrt to own shock
imp_y_u_mat = NaN(n_series,hor);
coef_y_u.Q = [1, zeros(1,n_uarlag-1)];
for is = 1:n_series
    coef_y_u.M = [uar_coef_mat(is,:); eye(n_uarlag-1),zeros(n_uarlag-1,1)];
    coef_y_u.G = [uar_ser_mat(is); zeros(n_uarlag-1,1)];
    imp_y_u_mat(is,:) = varirf_ss(coef_y_u,0,hor);   % irf to own shocks
end

% Transform IRFs so they are units of levels of series
for is = 1:n_series
    tc = tcodevec(is);
    tmp = units_to_levels(imp_y_u_mat(is,:)', tc);
    imp_y_u_mat(is,:) = tmp';
    for i_shock = 1:n_shocks
        tmp = units_to_levels(reshape(imp_y_fac_mat(is,i_shock,:),hor,1), tc);
        imp_y_fac_mat(is,i_shock,:) = tmp';
    end
end

% Compute Variance Components for each Shock
tmp = imp_y_u_mat.^2;
vcomp_y_u_mat = cumsum(tmp,2);
vcomp_y_fac_mat = imp_y_fac_mat.^2;

vtotal_y_fac_mat = NaN(n_series,hor);
for is = 1:n_series
    tmp = reshape(vcomp_y_fac_mat(is,:,:),n_shocks,hor);
    tmp = cumsum(tmp,2);
    vcomp_y_fac_mat(is,:,:) = tmp;
	vtotal_y_fac_mat(is,:) = sum(tmp,1);
end

vtotal_y_mat = vtotal_y_fac_mat + vcomp_y_u_mat;   % total variance by forecast horizon
vfrac_y_fac_mat = vtotal_y_fac_mat./vtotal_y_mat;
vfrac_y_comp_mat = NaN(n_series,n_shocks,hor);
for i = 1:n_shocks
    tmp = reshape(vcomp_y_fac_mat(:,i,:),n_series,hor);
    tmp = tmp./vtotal_y_mat;
    vfrac_y_comp_mat(:,i,:) = tmp;
end


% SAVE OUTPUT
out.imp_fac_mat     = imp_fac_mat;          % irf from factor VAR 
out.vcomp_fac_mat   = vcomp_fac_mat;        % variance component from factor VAR
out.vtotal_fac_mat  = vtotal_fac_mat;       % total variance from factor VAR
out.vfrac_fac_comp_mat = vfrac_fac_comp_mat;% fraction from factor VAR
out.imp_y_fac_mat     = imp_y_fac_mat;      % irf to factor shocks
out.imp_y_u_mat       = imp_y_u_mat ;       % irf for own shocks
out.vcomp_y_fac_mat   = vcomp_y_fac_mat;    % variance component from factors
out.vcomp_y_u_mat     = vcomp_y_u_mat;      % variance component from own shocks
out.vtotal_y_fac_mat  = vtotal_y_fac_mat;   % total variance from factors
out.vfrac_y_fac_mat   = vfrac_y_fac_mat;    % fraction from factor shocks
out.vfrac_y_comp_mat = vfrac_y_comp_mat;    % fraction 
out.vtotal_y_mat      = vtotal_y_mat;       % total variance
end