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