Description of bvar

0001 function bvar_forecast(nlags)
0002 % function bvar_forecast(nlags)
0003 % builds forecats for a bvar model
0004 %
0005 % INPUTS
0006 %    nlags:     number of lags for the bvar
0007 %
0008 % OUTPUTS
0009 %    none
0010 %
0011 % SPECIAL REQUIREMENTS
0012 %    none
0013 
0014 % Copyright (C) 2007-2010,2012 Dynare Team
0015 %
0016 % This file is part of Dynare.
0017 %
0018 % Dynare is free software: you can redistribute it and/or modify
0019 % it under the terms of the GNU General Public License as published by
0020 % the Free Software Foundation, either version 3 of the License, or
0021 % (at your option) any later version.
0022 %
0023 % Dynare is distributed in the hope that it will be useful,
0024 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0025 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0026 % GNU General Public License for more details.
0027 %
0028 % You should have received a copy of the GNU General Public License
0029 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0030 
0031 global options_ oo_ M_
0032 
0033 if options_.forecast == 0
0034     error('bvar_forecast: you must specify "forecast" option')
0035 end
0036 [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
0037 
0038 sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic);
0039 sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);
0040 
0041 S_inv_upper_chol = chol(inv(posterior.S));
0042 
0043 % Option 'lower' of chol() not available in old versions of
0044 % Matlab, so using transpose
0045 XXi_lower_chol = chol(posterior.XXi)';
0046 
0047 k = ny*nlags+nx;
0048 
0049 % Declaration of the companion matrix
0050 Companion_matrix = diag(ones(ny*(nlags-1),1),-ny);
0051 
0052 % Number of explosive VAR models sampled
0053 p = 0;
0054 % Loop counter initialization
0055 d = 0;
0056 while d <= options_.bvar_replic
0057     
0058     Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
0059 
0060     % Option 'lower' of chol() not available in old versions of
0061     % Matlab, so using transpose
0062     Sigma_lower_chol = chol(Sigma)';
0063 
0064     Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
0065     
0066     % All the eigenvalues of the companion matrix have to be on or inside the unit circle
0067     Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)'; 
0068     test = (abs(eig(Companion_matrix)));
0069     if any(test>1.0000000000001)
0070         p = p+1;
0071     end
0072     d = d+1;
0073     
0074     % Without shocks
0075     lags_data = forecast_data.initval;
0076     for t = 1:options_.forecast
0077         X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
0078         y = X * Phi;
0079         lags_data(1:end-1,:) = lags_data(2:end, :);
0080         lags_data(end,:) = y;
0081         sims_no_shock(t, :, d) = y;
0082     end
0083     
0084     % With shocks
0085     lags_data = forecast_data.initval;
0086     for t = 1:options_.forecast
0087         X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
0088         shock = (Sigma_lower_chol * randn(ny, 1))';
0089         y = X * Phi + shock;
0090         lags_data(1:end-1,:) = lags_data(2:end, :);
0091         lags_data(end,:) = y;
0092         sims_with_shocks(t, :, d) = y;
0093     end
0094 end
0095 
0096 if p > 0
0097     disp(' ')
0098     disp(['Some of the VAR models sampled from the posterior distribution'])
0099     disp(['were found to be explosive (' num2str(p/options_.bvar_replic) ' percent).'])
0100     disp(' ')
0101 end
0102 
0103 % Plot graphs
0104 sims_no_shock_mean = mean(sims_no_shock, 3);
0105 
0106 sort_idx = round((0.5 + [-options_.conf_sig, options_.conf_sig, 0]/2) * options_.bvar_replic);
0107 
0108 sims_no_shock_sort = sort(sims_no_shock, 3);
0109 sims_no_shock_down_conf = sims_no_shock_sort(:, :, sort_idx(1));
0110 sims_no_shock_up_conf = sims_no_shock_sort(:, :, sort_idx(2));
0111 sims_no_shock_median = sims_no_shock_sort(:, :, sort_idx(3));
0112 
0113 sims_with_shocks_sort = sort(sims_with_shocks, 3);
0114 sims_with_shocks_down_conf = sims_with_shocks_sort(:, :, sort_idx(1));
0115 sims_with_shocks_up_conf = sims_with_shocks_sort(:, :, sort_idx(2));
0116 
0117 dynare_graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});
0118 
0119 for i = 1:ny
0120     dynare_graph([ sims_no_shock_median(:, i) ...
0121                    sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
0122                    sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
0123                  options_.varobs(i, :));
0124 end
0125 
0126 dynare_graph_close;
0127 
0128 
0129 % Compute RMSE
0130 
0131 if ~isempty(forecast_data.realized_val)
0132     
0133     sq_err_cumul = zeros(1, ny);
0134     
0135     lags_data = forecast_data.initval;
0136     for t = 1:size(forecast_data.realized_val, 1)
0137         X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ];
0138         y = X * posterior.PhiHat;
0139         lags_data(1:end-1,:) = lags_data(2:end, :);
0140         lags_data(end,:) = y;
0141         sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2;
0142     end
0143     
0144     rmse = sqrt(sq_err_cumul / size(forecast_data.realized_val, 1));
0145     
0146     fprintf('RMSE of BVAR(%d):\n', nlags);
0147     
0148     for i = 1:size(options_.varobs, 1)
0149         fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i));
0150     end 
0151 end
0152 
0153 % Store results
0154 
0155 DirectoryName = [ M_.fname '/bvar_forecast' ];
0156 if ~isdir(DirectoryName)
0157     if ~isdir(M_.fname)
0158         mkdir(M_.fname);
0159     end
0160     mkdir(DirectoryName);
0161 end
0162 save([ DirectoryName '/simulations.mat'], 'sims_no_shock', 'sims_with_shocks');
0163 
0164 for i = 1:size(options_.varobs, 1)
0165     name = options_.varobs(i, :);
0166 
0167     sims = squeeze(sims_with_shocks(:,i,:));
0168     eval(['oo_.bvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']);
0169     eval(['oo_.bvar.forecast.with_shocks.Median.' name ' = median(sims, 2);']);
0170     eval(['oo_.bvar.forecast.with_shocks.Var.' name ' = var(sims, 0, 2);']);
0171     eval(['oo_.bvar.forecast.with_shocks.HPDsup.' name ' = sims_with_shocks_up_conf(:,i);']);
0172     eval(['oo_.bvar.forecast.with_shocks.HPDinf.' name ' = sims_with_shocks_down_conf(:,i);']);
0173 
0174     sims = squeeze(sims_no_shock(:,i,:));
0175     eval(['oo_.bvar.forecast.no_shock.Mean.' name ' = sims_no_shock_mean(:,i);']);
0176     eval(['oo_.bvar.forecast.no_shock.Median.' name ' = sims_no_shock_median(:,i);']);
0177     eval(['oo_.bvar.forecast.no_shock.Var.' name ' = var(sims, 0, 2);']);
0178     eval(['oo_.bvar.forecast.no_shock.HPDsup.' name ' = sims_no_shock_up_conf(:,i);']);
0179     eval(['oo_.bvar.forecast.no_shock.HPDinf.' name ' = sims_no_shock_down_conf(:,i);']);
0180 
0181     if exist('rmse')
0182         eval(['oo_.bvar.forecast.rmse.' name ' = rmse(i);']);
0183     end
0184 end
bvar_forecast

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
