bvar_irf.m 5.22 KB
Newer Older
1
function bvar_irf(nlags,identification)
2 3 4 5
% builds IRFs for a bvar model
%
% INPUTS
%    nlags            [integer]     number of lags for the bvar
6
%    identification   [string]      identification scheme ('Cholesky' or 'SquareRoot')
7 8 9 10 11 12 13
%
% OUTPUTS
%    none
%
% SPECIAL REQUIREMENTS
%    none

14
% Copyright (C) 2007-2012 Dynare Team
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
%
% 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/>.

31
global options_ oo_ M_
32

33 34 35
if nargin==1
    identification = 'Cholesky';
end
36

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
[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
57
    
58 59 60 61
    % 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);
62

63 64
    % Get the Autoregressive matrices from a matrix variate distribution.
    Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
65
    
66 67
    % Form the companion matrix.
    Companion_matrix(1:ny,:) = transpose(Phi(1:ny*nlags,:)); 
68
    
69 70 71 72 73
    % All the eigenvalues of the companion matrix have to be on or
    % inside the unit circle to rule out explosive time series.
    test = (abs(eig(Companion_matrix)));
    if any(test>1.0000000000001)
        p = p+1;
74 75
    end

76 77 78 79 80
    if strcmpi(identification,'Cholesky')
        StructuralMat = Sigma_lower_chol;
    elseif strcmpi(identification,'SquareRoot')
        StructuralMat = sqrtm(Sigma);
    end
81
    
82 83 84 85 86 87 88 89
    % Build the IRFs...
    lags_data = zeros(ny,ny*nlags) ;
    sampled_irfs(:,:,1,draw) = Sigma_lower_chol ;
    lags_data(:,1:ny) = Sigma_lower_chol ;
    for t=2:options_.irf
        sampled_irfs(:,:,t,draw) = lags_data(:,:)*Phi(1:ny*nlags,:) ;
        lags_data(:,ny+1:end) = lags_data(:,1:end-ny) ;
        lags_data(:,1:ny) = sampled_irfs(:,:,t,draw) ;
90 91
    end
    
92 93 94
end

if p > 0
95
    skipline()
96 97
    disp(['Some of the VAR models sampled from the posterior distribution'])
    disp(['were found to be explosive (' int2str(p) ' samples).'])
98
    skipline()
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
end

posterior_mean_irfs = mean(sampled_irfs,4);
posterior_variance_irfs = var(sampled_irfs, 1, 4);

sorted_irfs = sort(sampled_irfs,4);
sort_idx = round((0.5 + [-options_.conf_sig, options_.conf_sig, .0]/2) * options_.bvar_replic);

posterior_down_conf_irfs = sorted_irfs(:,:,:,sort_idx(1));
posterior_up_conf_irfs = sorted_irfs(:,:,:,sort_idx(2));
posterior_median_irfs = sorted_irfs(:,:,:,sort_idx(3));   

number_of_columns = fix(sqrt(ny));
number_of_rows = ceil(ny / number_of_columns) ;

% Plots of the IRFs
for shock=1:ny
    figure('Name',['Posterior BVAR Impulse Responses (shock in equation ' int2str(shock) ').']);
    for variable=1:ny
        subplot(number_of_rows,number_of_columns,variable);
119
        h1 = area(1:options_.irf,squeeze(posterior_up_conf_irfs(shock,variable,:)),'FaceColor',[.9 .9 .9],'BaseValue',min([min(posterior_up_conf_irfs(shock,variable,:)),min(posterior_down_conf_irfs(shock,variable,:))]));
120
        hold on
121
        h2 = area(1:options_.irf,squeeze(posterior_down_conf_irfs(shock,variable,:)),'FaceColor',[1 1 1],'BaseValue',min([min(posterior_up_conf_irfs(shock,variable,:)),min(posterior_down_conf_irfs(shock,variable,:))]));
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
        plot(1:options_.irf,squeeze(posterior_median_irfs(shock,variable,:)),'-k','linewidth',2)
        axis tight
        hold off
    end
end

% Save intermediate results
DirectoryName = [ M_.fname '/bvar_irf' ];
if ~isdir(DirectoryName)
    mkdir('.',DirectoryName);
end
save([ DirectoryName '/simulations.mat'], 'sampled_irfs');

% Save results in oo_
for i=1:ny
137
    shock_name = options_.varobs{i};
138
    for j=1:ny
139
        variable_name = options_.varobs{j};
140 141 142 143 144
        eval(['oo_.bvar.irf.Mean.' variable_name '.' shock_name ' = posterior_mean_irfs(' int2str(j) ',' int2str(i) ',:);'])
        eval(['oo_.bvar.irf.Median.' variable_name '.' shock_name ' = posterior_median_irfs(' int2str(j) ',' int2str(i) ',:);'])
        eval(['oo_.bvar.irf.Var.' variable_name '.' shock_name ' = posterior_variance_irfs(' int2str(j) ',' int2str(i) ',:);'])
        eval(['oo_.bvar.irf.Upper_bound.' variable_name '.' shock_name ' = posterior_up_conf_irfs(' int2str(j) ',' int2str(i) ',:);'])
        eval(['oo_.bvar.irf.Lower_bound.' variable_name '.' shock_name ' = posterior_down_conf_irfs(' int2str(j) ',' int2str(i) ',:);'])
145
    end
146
end