1function bvar_forecast(nlags) 2% function bvar_forecast(nlags) 3% builds forecats for a bvar model 4% 5% INPUTS 6% nlags: number of lags for the bvar 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 options_.forecast == 0 34 error('bvar_forecast: you must specify "forecast" option') 35end 36[ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags); 37 38sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic); 39sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic); 40 41S_inv_upper_chol = chol(inv(posterior.S)); 42 43% Option 'lower' of chol() not available in old versions of 44% Matlab, so using transpose 45XXi_lower_chol = chol(posterior.XXi)'; 46 47k = ny*nlags+nx; 48 49% Declaration of the companion matrix 50Companion_matrix = diag(ones(ny*(nlags-1),1),-ny); 51 52% Number of explosive VAR models sampled 53p = 0; 54% Loop counter initialization 55d = 0; 56while d <= options_.bvar_replic 57 58 Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol); 59 60 % Option 'lower' of chol() not available in old versions of 61 % Matlab, so using transpose 62 Sigma_lower_chol = chol(Sigma)'; 63 64 Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol); 65 66 % All the eigenvalues of the companion matrix have to be on or inside the unit circle 67 Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)'; 68 test = (abs(eig(Companion_matrix))); 69 if any(test>1.0000000000001) 70 p = p+1; 71 end 72 d = d+1; 73 74 % Without shocks 75 lags_data = forecast_data.initval; 76 for t = 1:options_.forecast 77 X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ]; 78 y = X * Phi; 79 lags_data(1:end-1,:) = lags_data(2:end, :); 80 lags_data(end,:) = y; 81 sims_no_shock(t, :, d) = y; 82 end 83 84 % With shocks 85 lags_data = forecast_data.initval; 86 for t = 1:options_.forecast 87 X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ]; 88 shock = (Sigma_lower_chol * randn(ny, 1))'; 89 y = X * Phi + shock; 90 lags_data(1:end-1,:) = lags_data(2:end, :); 91 lags_data(end,:) = y; 92 sims_with_shocks(t, :, d) = y; 93 end 94end 95 96if p > 0 97 skipline() 98 disp(['Some of the VAR models sampled from the posterior distribution']) 99 disp(['were found to be explosive (' num2str(p/options_.bvar_replic) ' percent).']) 100 skipline() 101end 102 103% Plot graphs 104sims_no_shock_mean = mean(sims_no_shock, 3); 105 106sort_idx = round((0.5 + [-options_.bvar.conf_sig, options_.bvar.conf_sig, 0]/2) * options_.bvar_replic); 107 108sims_no_shock_sort = sort(sims_no_shock, 3); 109sims_no_shock_down_conf = sims_no_shock_sort(:, :, sort_idx(1)); 110sims_no_shock_up_conf = sims_no_shock_sort(:, :, sort_idx(2)); 111sims_no_shock_median = sims_no_shock_sort(:, :, sort_idx(3)); 112 113sims_with_shocks_sort = sort(sims_with_shocks, 3); 114sims_with_shocks_down_conf = sims_with_shocks_sort(:, :, sort_idx(1)); 115sims_with_shocks_up_conf = sims_with_shocks_sort(:, :, sort_idx(2)); 116 117OutputDirectoryName = CheckPath('graphs',M_.fname); 118 119dyn_graph=dynare_graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'}); 120 121for i = 1:ny 122 dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ... 123 sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ... 124 sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ... 125 options_.varobs{i}); 126end 127 128dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_BVAR_forecast_',num2str(nlags)],options_.nodisplay,options_.graph_format) 129 130% Compute RMSE 131 132if ~isempty(forecast_data.realized_val) 133 134 sq_err_cumul = zeros(1, ny); 135 136 lags_data = forecast_data.initval; 137 for t = 1:size(forecast_data.realized_val, 1) 138 X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ]; 139 y = X * posterior.PhiHat; 140 lags_data(1:end-1,:) = lags_data(2:end, :); 141 lags_data(end,:) = y; 142 sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2; 143 end 144 145 rmse = sqrt(sq_err_cumul / size(forecast_data.realized_val, 1)); 146 147 fprintf('RMSE of BVAR(%d):\n', nlags); 148 149 for i = 1:length(options_.varobs) 150 fprintf('%s: %10.4f\n', options_.varobs{i}, rmse(i)); 151 end 152end 153 154% Store results 155 156DirectoryName = [ M_.fname '/bvar_forecast' ]; 157if ~isdir(DirectoryName) 158 if ~isdir(M_.fname) 159 mkdir(M_.fname); 160 end 161 mkdir(DirectoryName); 162end 163save([ DirectoryName '/simulations.mat'], 'sims_no_shock', 'sims_with_shocks'); 164 165for i = 1:length(options_.varobs) 166 name = options_.varobs{i}; 167 168 sims = squeeze(sims_with_shocks(:,i,:)); 169 eval(['oo_.bvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']); 170 eval(['oo_.bvar.forecast.with_shocks.Median.' name ' = median(sims, 2);']); 171 eval(['oo_.bvar.forecast.with_shocks.Var.' name ' = var(sims, 0, 2);']); 172 eval(['oo_.bvar.forecast.with_shocks.HPDsup.' name ' = sims_with_shocks_up_conf(:,i);']); 173 eval(['oo_.bvar.forecast.with_shocks.HPDinf.' name ' = sims_with_shocks_down_conf(:,i);']); 174 175 sims = squeeze(sims_no_shock(:,i,:)); 176 eval(['oo_.bvar.forecast.no_shock.Mean.' name ' = sims_no_shock_mean(:,i);']); 177 eval(['oo_.bvar.forecast.no_shock.Median.' name ' = sims_no_shock_median(:,i);']); 178 eval(['oo_.bvar.forecast.no_shock.Var.' name ' = var(sims, 0, 2);']); 179 eval(['oo_.bvar.forecast.no_shock.HPDsup.' name ' = sims_no_shock_up_conf(:,i);']); 180 eval(['oo_.bvar.forecast.no_shock.HPDinf.' name ' = sims_no_shock_down_conf(:,i);']); 181 182 if exist('rmse') 183 eval(['oo_.bvar.forecast.rmse.' name ' = rmse(i);']); 184 end 185end 186