From eb0ccd8f5b514bc5fb7854fddd4bd19a462fe69c Mon Sep 17 00:00:00 2001 From: adjemian Date: Fri, 12 Sep 2008 10:14:36 +0000 Subject: [PATCH] v4.1: Added posterior IRFs for bvar a la Sims (not yet ready. Results are not saved in oo_, more general identification conditions should be allowed, and the code has to be tested). git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2059 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/bvar_irf.m | 121 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100644 matlab/bvar_irf.m diff --git a/matlab/bvar_irf.m b/matlab/bvar_irf.m new file mode 100644 index 000000000..eaba2cdb2 --- /dev/null +++ b/matlab/bvar_irf.m @@ -0,0 +1,121 @@ +function bvar_irf(nlags,varagin) +% builds IRFs for a bvar model +% +% INPUTS +% nlags [integer] number of lags for the bvar +% identification [string] identification scheme ('cholesky') +% +% OUTPUTS +% none +% +% SPECIAL REQUIREMENTS +% none + +% Copyright (C) 2007-2008 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 . + + global options_ oo_ M_ + + options_ = set_default_option(options_, 'bvar_replic', 2000); + + [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. + test = (abs(eig(Companion_matrix))); + if any(test>1.0000000000001) + p = p+1; + end + + % 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) ; + end + + end + + if p > 0 + disp(' ') + disp(['Some of the VAR models sampled from the posterior distribution']) + disp(['were found to be explosive (' int2str(p) ' samples).']) + disp(' ') + 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) ; + + 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); + h1 = area(1:options_.irf,squeeze(posterior_up_conf_irfs(shock,variable,:))); + set(h1,'BaseValue',min([min(posterior_up_conf_irfs(shock,variable,:)),min(posterior_down_conf_irfs(shock,variable,:))])) + set(h1,'FaceColor',[.9 .9 .9]) + hold on + h2 = area(1:options_.irf,squeeze(posterior_down_conf_irfs(shock,variable,:))); + set(h2,'BaseValue',min([min(posterior_up_conf_irfs(shock,variable,:)),min(posterior_down_conf_irfs(shock,variable,:))])) + set(h2,'FaceColor',[1 1 1]) + plot(1:options_.irf,squeeze(posterior_median_irfs(shock,variable,:)),'-k','linewidth',2) + axis tight + hold off + end + end \ No newline at end of file