bvar_irf.m

function bvar_irf(nlags,identification)
% builds IRFs for a bvar model
%
% INPUTS
%    nlags            [integer]     number of lags for the bvar
%    identification   [string]      identification scheme ('Cholesky' or 'SquareRoot')
%
% OUTPUTS
%    none
%
% SPECIAL REQUIREMENTS
%    none

% Copyright (C) 2007-2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

global options_ oo_ M_

if nargin==1
    identification = 'Cholesky';
end

[ny, nx, posterior, prior] = bvar_toolbox(nlags);

S_inv_upper_chol = chol(inv(posterior.S));

% Option 'lower' of chol() not available in old versions of
% Matlab, so using transpose
XXi_lower_chol = chol(posterior.XXi)';

k = ny*nlags+nx;

% Declaration of the companion matrix
Companion_matrix = diag(ones(ny*(nlags-1),1),-ny);

% Number of explosive VAR models sampled
p = 0;

% Initialize a four dimensional array.
sampled_irfs = NaN(ny, ny, options_.irf, options_.bvar_replic);

for draw=1:options_.bvar_replic
    
    % Get a covariance matrix from an inverted Wishart distribution.
    Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
    Sigma_upper_chol = chol(Sigma);
    Sigma_lower_chol = transpose(Sigma_upper_chol);

    % Get the Autoregressive matrices from a matrix variate distribution.
    Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
    
    % Form the companion matrix.
    Companion_matrix(1:ny,:) = transpose(Phi(1:ny*nlags,:)); 
    
    % All the eigenvalues of the companion matrix have to be on or
    % inside the unit circle to rule out explosive time series.