forcst_unc

 function [mean,intval1,intval2]=forcst_unc(y0,var_list)
 computes forecasts with parameter uncertainty

 INPUTS
   y0: matrix of initial values
   var_list: list of variables to be forecasted

 OUTPUTS
   none

 ALGORITHM
   uses antithetic draws for the shocks

 SPECIAL REQUIREMENTS
   None.

0001 function forcst_unc(y0,var_list)
0002 % function [mean,intval1,intval2]=forcst_unc(y0,var_list)
0003 % computes forecasts with parameter uncertainty
0004 %
0005 % INPUTS
0006 %   y0: matrix of initial values
0007 %   var_list: list of variables to be forecasted
0008 %
0009 % OUTPUTS
0010 %   none
0011 %
0012 % ALGORITHM
0013 %   uses antithetic draws for the shocks
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %   None.
0017 
0018 % Copyright (C) 2006-2011 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 global M_ options_ oo_ estim_params_ bayestopt_
0036 
0037 % setting up estim_params_
0038 [xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_,M_);
0039 
0040 options_.TeX = 0;
0041 options_.nograph = 0;
0042 plot_priors(bayestopt_,M_,options_);
0043 
0044 % workspace initialization
0045 if isempty(var_list)
0046     var_list = M_.endo_names(1:M_.orig_endo_nbr,:);
0047 end
0048 n = size(var_list,1);
0049 
0050 periods = options_.forecast;
0051 exo_nbr = M_.exo_nbr;
0052 replic = options_.replic;
0053 order = options_.order;
0054 maximum_lag = M_.maximum_lag;
0055 %  params = prior_draw(1);
0056 params = rndprior(bayestopt_);
0057 set_parameters(params);
0058 % eliminate shocks with 0 variance
0059 i_exo_var = setdiff([1:exo_nbr],find(diag(M_.Sigma_e) == 0));
0060 nx = length(i_exo_var);
0061 
0062 ex0 = zeros(periods,exo_nbr);
0063 yf1 = zeros(periods+M_.maximum_lag,n,replic);
0064 
0065 % loops on parameter values
0066 m1 = 0;
0067 m2 = 0;
0068 for i=1:replic
0069     % draw parameter values from the prior
0070     % params = prior_draw(0);
0071     params = rndprior(bayestopt_);
0072     set_parameters(params);
0073     % solve the model
0074     [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0075     % discard problematic cases
0076     if info
0077         continue
0078     end
0079     % compute forecast with zero shocks
0080     m1 = m1+1;
0081     yf1(:,:,m1) = simult_(y0,dr,ex0,order)';
0082     % compute forecast with antithetic shocks
0083     chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var));
0084     ex(:,i_exo_var) = randn(periods,nx)*chol_S;
0085     m2 = m2+1;
0086     yf2(:,:,m2) = simult_(y0,dr,ex,order)';
0087     m2 = m2+1;
0088     yf2(:,:,m2) = simult_(y0,dr,-ex,order)';
0089 end
0090 
0091 oo_.forecast.accept_rate = (replic-m1)/replic;
0092 
0093 if options_.noprint == 0 && m1 < replic
0094     disp(' ')
0095     disp(' ')
0096     disp('FORECASTING WITH PARAMETER UNCERTAINTY:')
0097     disp(sprintf(['The model  couldn''t be solved for %f%% of the parameter' ...
0098                   ' values'],100*oo_.forecast.accept_rate))
0099     disp(' ')
0100     disp(' ')
0101 end
0102 
0103 % compute moments
0104 yf1 = yf1(:,:,1:m1);
0105 yf2 = yf2(:,:,1:m2);
0106 
0107 yf_mean = mean(yf1,3);
0108 
0109 yf1 = sort(yf1,3);
0110 yf2 = sort(yf2,3);
0111 
0112 sig_lev = options_.conf_sig;
0113 k1 = round((0.5+[-sig_lev, sig_lev]/2)*replic);
0114 % make sure that lower bound is at least the first element
0115 if k1(2) == 0
0116     k1(2) = 1;
0117 end
0118 k2 = round((1+[-sig_lev, sig_lev])*replic);
0119 % make sure that lower bound is at least the first element
0120 if k2(2) == 0
0121     k2(2) = 1;
0122 end
0123 
0124 % compute shock uncertainty around forecast with mean prior
0125 set_parameters(bayestopt_.p1);
0126 [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0127 [yf3,yf3_intv] = forcst(dr,y0,periods,var_list);
0128 yf3_1 = yf3'-[zeros(maximum_lag,n); yf3_intv];
0129 yf3_2 = yf3'+[zeros(maximum_lag,n); yf3_intv];
0130 
0131 % graphs
0132 
0133 dynare_graph_init('Forecasts type I',n,{'b-' 'g-' 'g-' 'r-' 'r-'});
0134 for i=1:n
0135     dynare_graph([yf_mean(:,i) squeeze(yf1(:,i,k1)) squeeze(yf2(:,i,k2))],...
0136                  var_list(i,:));
0137 end
0138 dynare_graph_close;
0139 
0140 dynare_graph_init('Forecasts type II',n,{'b-' 'k-' 'k-' 'r-' 'r-'});
0141 for i=1:n
0142     dynare_graph([yf_mean(:,i) yf3_1(:,i) yf3_2(:,i) squeeze(yf2(:,i,k2))],...
0143                  var_list(i,:));
0144 end
0145 dynare_graph_close;
0146 
0147 
0148 % saving results
0149 save_results(yf_mean,'oo_.forecast.mean.',var_list);
0150 save_results(yf1(:,:,k1(1)),'oo_.forecast.HPDinf.',var_list);
0151 save_results(yf1(:,:,k1(2)),'oo_.forecast.HPDsup.',var_list);
0152 save_results(yf2(:,:,k2(1)),'oo_.forecast.HPDTotalinf.',var_list);
0153 save_results(yf2(:,:,k2(2)),'oo_.forecast.HPDTotalsup.',var_list);