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