2007-06-26 00:42:30 +02:00
|
|
|
function bvar_forecast(nlags)
|
2008-03-03 11:52:19 +01:00
|
|
|
% function bvar_forecast(nlags)
|
|
|
|
% builds forecats for a bvar model
|
|
|
|
%
|
|
|
|
% INPUTS
|
|
|
|
% nlags: number of lags for the bvar
|
|
|
|
%
|
|
|
|
% OUTPUTS
|
|
|
|
% none
|
|
|
|
%
|
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
% none
|
2008-08-01 14:40:33 +02:00
|
|
|
|
2011-02-04 17:27:33 +01:00
|
|
|
% Copyright (C) 2007-2010 Dynare Team
|
2008-08-01 14:40:33 +02:00
|
|
|
%
|
|
|
|
% 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 <http://www.gnu.org/licenses/>.
|
2008-03-03 11:52:19 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
global options_ oo_ M_
|
|
|
|
|
|
|
|
options_ = set_default_option(options_, 'bvar_replic', 2000);
|
|
|
|
if options_.forecast == 0
|
|
|
|
error('bvar_forecast: you must specify "forecast" option')
|
|
|
|
end
|
|
|
|
[ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
|
|
|
|
|
|
|
|
sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic);
|
|
|
|
sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);
|
|
|
|
|
|
|
|
S_inv_upper_chol = chol(inv(posterior.S));
|
2007-09-21 14:17:56 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
% 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;
|
|
|
|
% Loop counter initialization
|
|
|
|
d = 0;
|
|
|
|
while d <= options_.bvar_replic
|
2007-06-26 00:42:30 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
|
2007-07-05 14:49:26 +02:00
|
|
|
|
|
|
|
% Option 'lower' of chol() not available in old versions of
|
|
|
|
% Matlab, so using transpose
|
2009-12-16 18:17:34 +01:00
|
|
|
Sigma_lower_chol = chol(Sigma)';
|
|
|
|
|
|
|
|
Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
|
2007-06-26 00:42:30 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
% All the eigenvalues of the companion matrix have to be on or inside the unit circle
|
|
|
|
Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)';
|
|
|
|
test = (abs(eig(Companion_matrix)));
|
|
|
|
if any(test>1.0000000000001)
|
|
|
|
p = p+1;
|
2007-06-26 00:42:30 +02:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
d = d+1;
|
2007-09-21 14:17:56 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
% Without shocks
|
|
|
|
lags_data = forecast_data.initval;
|
|
|
|
for t = 1:options_.forecast
|
|
|
|
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
|
|
|
|
y = X * Phi;
|
|
|
|
lags_data(1:end-1,:) = lags_data(2:end, :);
|
|
|
|
lags_data(end,:) = y;
|
|
|
|
sims_no_shock(t, :, d) = y;
|
2007-12-07 15:27:37 +01:00
|
|
|
end
|
2007-09-21 14:17:56 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
% With shocks
|
|
|
|
lags_data = forecast_data.initval;
|
|
|
|
for t = 1:options_.forecast
|
|
|
|
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
|
|
|
|
shock = (Sigma_lower_chol * randn(ny, 1))';
|
|
|
|
y = X * Phi + shock;
|
|
|
|
lags_data(1:end-1,:) = lags_data(2:end, :);
|
|
|
|
lags_data(end,:) = y;
|
|
|
|
sims_with_shocks(t, :, d) = y;
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
if p > 0
|
|
|
|
disp(' ')
|
|
|
|
disp(['Some of the VAR models sampled from the posterior distribution'])
|
|
|
|
disp(['were found to be explosive (' num2str(p/options_.bvar_replic) ' percent).'])
|
|
|
|
disp(' ')
|
|
|
|
end
|
2007-06-26 00:42:30 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
% Plot graphs
|
|
|
|
sims_no_shock_mean = mean(sims_no_shock, 3);
|
|
|
|
|
|
|
|
sort_idx = round((0.5 + [-options_.conf_sig, options_.conf_sig, 0]/2) * options_.bvar_replic);
|
|
|
|
|
|
|
|
sims_no_shock_sort = sort(sims_no_shock, 3);
|
|
|
|
sims_no_shock_down_conf = sims_no_shock_sort(:, :, sort_idx(1));
|
|
|
|
sims_no_shock_up_conf = sims_no_shock_sort(:, :, sort_idx(2));
|
|
|
|
sims_no_shock_median = sims_no_shock_sort(:, :, sort_idx(3));
|
|
|
|
|
|
|
|
sims_with_shocks_sort = sort(sims_with_shocks, 3);
|
|
|
|
sims_with_shocks_down_conf = sims_with_shocks_sort(:, :, sort_idx(1));
|
|
|
|
sims_with_shocks_up_conf = sims_with_shocks_sort(:, :, sort_idx(2));
|
|
|
|
|
|
|
|
dynare_graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});
|
|
|
|
|
|
|
|
for i = 1:ny
|
|
|
|
dynare_graph([ sims_no_shock_median(:, i) ...
|
|
|
|
sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
|
|
|
|
sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
|
|
|
|
options_.varobs(i, :));
|
|
|
|
end
|
|
|
|
|
|
|
|
dynare_graph_close;
|
|
|
|
|
|
|
|
|
|
|
|
% Compute RMSE
|
|
|
|
|
|
|
|
if ~isempty(forecast_data.realized_val)
|
2007-06-26 00:42:30 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
sq_err_cumul = zeros(1, ny);
|
|
|
|
|
|
|
|
lags_data = forecast_data.initval;
|
|
|
|
for t = 1:size(forecast_data.realized_val, 1)
|
|
|
|
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ];
|
|
|
|
y = X * posterior.PhiHat;
|
|
|
|
lags_data(1:end-1,:) = lags_data(2:end, :);
|
|
|
|
lags_data(end,:) = y;
|
|
|
|
sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2;
|
2007-06-26 00:42:30 +02:00
|
|
|
end
|
|
|
|
|
2010-01-08 10:31:57 +01:00
|
|
|
rmse = sqrt(sq_err_cumul / size(forecast_data.realized_val, 1));
|
2007-06-26 00:42:30 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
fprintf('RMSE of BVAR(%d):\n', nlags);
|
2007-06-26 00:42:30 +02:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
for i = 1:size(options_.varobs, 1)
|
|
|
|
fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i));
|
|
|
|
end
|
|
|
|
end
|
2007-12-07 15:27:37 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
% Store results
|
2007-12-07 15:27:37 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
DirectoryName = [ M_.fname '/bvar_forecast' ];
|
|
|
|
if ~isdir(DirectoryName)
|
|
|
|
if ~isdir(M_.fname)
|
|
|
|
mkdir(M_.fname);
|
2007-12-07 15:27:37 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
mkdir(DirectoryName);
|
|
|
|
end
|
|
|
|
save([ DirectoryName '/simulations.mat'], 'sims_no_shock', 'sims_with_shocks');
|
2007-12-07 15:27:37 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
for i = 1:size(options_.varobs, 1)
|
|
|
|
name = options_.varobs(i, :);
|
|
|
|
|
|
|
|
sims = squeeze(sims_with_shocks(:,i,:));
|
|
|
|
eval(['oo_.bvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']);
|
|
|
|
eval(['oo_.bvar.forecast.with_shocks.Median.' name ' = median(sims, 2);']);
|
|
|
|
eval(['oo_.bvar.forecast.with_shocks.Var.' name ' = var(sims, 0, 2);']);
|
|
|
|
eval(['oo_.bvar.forecast.with_shocks.HPDsup.' name ' = sims_with_shocks_up_conf(:,i);']);
|
|
|
|
eval(['oo_.bvar.forecast.with_shocks.HPDinf.' name ' = sims_with_shocks_down_conf(:,i);']);
|
|
|
|
|
|
|
|
sims = squeeze(sims_no_shock(:,i,:));
|
|
|
|
eval(['oo_.bvar.forecast.no_shock.Mean.' name ' = sims_no_shock_mean(:,i);']);
|
|
|
|
eval(['oo_.bvar.forecast.no_shock.Median.' name ' = sims_no_shock_median(:,i);']);
|
|
|
|
eval(['oo_.bvar.forecast.no_shock.Var.' name ' = var(sims, 0, 2);']);
|
|
|
|
eval(['oo_.bvar.forecast.no_shock.HPDsup.' name ' = sims_no_shock_up_conf(:,i);']);
|
|
|
|
eval(['oo_.bvar.forecast.no_shock.HPDinf.' name ' = sims_no_shock_down_conf(:,i);']);
|
|
|
|
|
|
|
|
if exist('rmse')
|
|
|
|
eval(['oo_.bvar.forecast.rmse.' name ' = rmse(i);']);
|
2007-12-07 15:27:37 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|