1function bvar_irf(nlags,identification) 2% builds IRFs for a bvar model 3% 4% INPUTS 5% nlags [integer] number of lags for the bvar 6% identification [string] identification scheme ('Cholesky' or 'SquareRoot') 7% 8% OUTPUTS 9% none 10% 11% SPECIAL REQUIREMENTS 12% none 13 14% Copyright (C) 2007-2017 Dynare Team 15% 16% This file is part of Dynare. 17% 18% Dynare is free software: you can redistribute it and/or modify 19% it under the terms of the GNU General Public License as published by 20% the Free Software Foundation, either version 3 of the License, or 21% (at your option) any later version. 22% 23% Dynare is distributed in the hope that it will be useful, 24% but WITHOUT ANY WARRANTY; without even the implied warranty of 25% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 26% GNU General Public License for more details. 27% 28% You should have received a copy of the GNU General Public License 29% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 30 31global options_ oo_ M_ 32 33if nargin==1 34 identification = 'Cholesky'; 35end 36 37[ny, nx, posterior, prior] = bvar_toolbox(nlags); 38 39S_inv_upper_chol = chol(inv(posterior.S)); 40 41% Option 'lower' of chol() not available in old versions of 42% Matlab, so using transpose 43XXi_lower_chol = chol(posterior.XXi)'; 44 45k = ny*nlags+nx; 46 47% Declaration of the companion matrix 48Companion_matrix = diag(ones(ny*(nlags-1),1),-ny); 49 50% Number of explosive VAR models sampled 51p = 0; 52 53% Initialize a four dimensional array. 54sampled_irfs = NaN(ny, ny, options_.irf, options_.bvar_replic); 55 56for draw=1:options_.bvar_replic 57 58 % Get a covariance matrix from an inverted Wishart distribution. 59 Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol); 60 Sigma_upper_chol = chol(Sigma); 61 Sigma_lower_chol = transpose(Sigma_upper_chol); 62 63 % Get the Autoregressive matrices from a matrix variate distribution. 64 Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol); 65 66 % Form the companion matrix. 67 Companion_matrix(1:ny,:) = transpose(Phi(1:ny*nlags,:)); 68 69 % All the eigenvalues of the companion matrix have to be on or 70 % inside the unit circle to rule out explosive time series. 71 test = (abs(eig(Companion_matrix))); 72 if any(test>1.0000000000001) 73 p = p+1; 74 end 75 76 if strcmpi(identification,'Cholesky') 77 StructuralMat = Sigma_lower_chol; 78 elseif strcmpi(identification,'SquareRoot') 79 StructuralMat = sqrtm(Sigma); 80 end 81 82 % Build the IRFs... 83 lags_data = zeros(ny,ny*nlags) ; 84 sampled_irfs(:,:,1,draw) = Sigma_lower_chol ; 85 lags_data(:,1:ny) = Sigma_lower_chol ; 86 for t=2:options_.irf 87 sampled_irfs(:,:,t,draw) = lags_data(:,:)*Phi(1:ny*nlags,:) ; 88 lags_data(:,ny+1:end) = lags_data(:,1:end-ny) ; 89 lags_data(:,1:ny) = sampled_irfs(:,:,t,draw) ; 90 end 91 92end 93 94if p > 0 95 skipline() 96 disp(['Some of the VAR models sampled from the posterior distribution']) 97 disp(['were found to be explosive (' int2str(p) ' samples).']) 98 skipline() 99end 100 101posterior_mean_irfs = mean(sampled_irfs,4); 102posterior_variance_irfs = var(sampled_irfs, 1, 4); 103 104sorted_irfs = sort(sampled_irfs,4); 105sort_idx = round((0.5 + [-options_.bvar.conf_sig, options_.bvar.conf_sig, .0]/2) * options_.bvar_replic); 106 107posterior_down_conf_irfs = sorted_irfs(:,:,:,sort_idx(1)); 108posterior_up_conf_irfs = sorted_irfs(:,:,:,sort_idx(2)); 109posterior_median_irfs = sorted_irfs(:,:,:,sort_idx(3)); 110 111number_of_columns = fix(sqrt(ny)); 112number_of_rows = ceil(ny / number_of_columns) ; 113 114% Plots of the IRFs 115for shock=1:ny 116 figure('Name',['Posterior BVAR Impulse Responses (shock in equation ' int2str(shock) ').']); 117 for variable=1:ny 118 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 plot(1:options_.irf,squeeze(posterior_median_irfs(shock,variable,:)),'-k','linewidth',2) 123 axis tight 124 hold off 125 end 126end 127 128% Save intermediate results 129DirectoryName = [ M_.fname '/bvar_irf' ]; 130if ~isdir(DirectoryName) 131 mkdir('.',DirectoryName); 132end 133save([ DirectoryName '/simulations.mat'], 'sampled_irfs'); 134 135% Save results in oo_ 136for i=1:ny 137 shock_name = options_.varobs{i}; 138 for j=1:ny 139 variable_name = options_.varobs{j}; 140 eval(['oo_.bvar.irf.Mean.' variable_name '.' shock_name ' = posterior_mean_irfs(' int2str(j) ',' int2str(i) ',:);']) 141 eval(['oo_.bvar.irf.Median.' variable_name '.' shock_name ' = posterior_median_irfs(' int2str(j) ',' int2str(i) ',:);']) 142 eval(['oo_.bvar.irf.Var.' variable_name '.' shock_name ' = posterior_variance_irfs(' int2str(j) ',' int2str(i) ',:);']) 143 eval(['oo_.bvar.irf.Upper_bound.' variable_name '.' shock_name ' = posterior_up_conf_irfs(' int2str(j) ',' int2str(i) ',:);']) 144 eval(['oo_.bvar.irf.Lower_bound.' variable_name '.' shock_name ' = posterior_down_conf_irfs(' int2str(j) ',' int2str(i) ',:);']) 145 end 146end