function vardecomp=get_var_decomp(L_companion,Pi_companion,Steps,n_lags,n_vars,variable_names)
% function vardecomp=get_var_decomp(L_companion,Pi_companion,Steps,n_lags,n_vars,variable_names)
% Inputs:
%   L_companion     [matrix]    Lower Cholesky of covariance matrix of
%                               companion form
%   Pi_companion    [matrix]    companion form VAR matrices
%   Steps           [vector]    horizon for variance decomposition
%   n_lags          [scalar]    number of lags
%   n_vars          [scalar]    number of variables in VAR
%   variable_names  [char]      array of variable names
% Output:
%   vardecomp       [structure] if variable_names is provided
%   vardecomp       [matrix]    if variable_names is not provided
%
% Notes: variables along second dimension; shocks along third dimension

%
% Copyright (C) 2015-2016 Benjamin Born and Johannes Pfeifer
%
%  This is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
%
%  It is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  For a copy of the GNU General Public License,
%  see <http://www.gnu.org/licenses/>.


n_horizon = length(Steps);
vardec_array = NaN(n_horizon,n_vars,n_vars);
%% do variance decomposition.
variance = zeros(n_horizon,n_vars);   % Temporary array to store variances;
for shock_iter = 1:n_vars   % Iterate over orthogonalized shocks:
    % Matrix that identifies the shock:
    W_selector = zeros(n_vars*n_lags);
    W_selector(shock_iter,shock_iter)=1;
    % Initialize variance of forecast error in t=1
    v_1 = L_companion*W_selector*L_companion';
    var_i = v_1;
    % Store just the variance of vari:
    vardec_array(1,:,shock_iter) =  diag(var_i(1:n_vars,1:n_vars))';
    
    % Compute recursively variance of t-step-ahead forecast error
    for t = 2:n_horizon
        var_i = v_1 + Pi_companion*var_i*Pi_companion';                                              % Update covariance matrix of forecast errors
        vardec_array(t,:,shock_iter) =  diag(var_i(1:n_vars,1:n_vars))';
    end
    % Update overall variance of t-step-ahead forecast error:
    variance = variance + squeeze(vardec_array(:,:,shock_iter));
end

% Now put this information as a percentage of the forecast error variance
for var_iter=1:n_vars        % For each series
    for shock_iter=1:n_vars    % For each shock
        vardec_array(:,var_iter,shock_iter) =  100*vardec_array(:,var_iter,shock_iter)./variance(:,var_iter);
    end
end

if max(max(abs(sum(vardec_array,3)-100)))>1e-6
    error('Variance decomposition does not add up')
end
if nargin>5
    vardecomp.vardec_array=vardec_array;
    vardecomp.variable_names=variable_names;
else
vardecomp=vardec_array;
end
end
