1function imcforecast(constrained_paths, constrained_vars, options_cond_fcst) 2 3% Computes conditional forecasts. 4% 5% INPUTS 6% - consnstrained_paths [double] m*p array, where m is the number of constrained endogenous variables and p is the number of constrained periods. 7% - constrained_vars [integer] m*1 array, indices in M_.endo_names of the constrained variables. 8% - options_cond_fcst [structure] containing the options. The fields are: 9% 10% + replic [integer] scalar, number of monte carlo simulations. 11% + parameter_set [char] values of the estimated parameters: 12% 'posterior_mode', 13% 'posterior_mean', 14% 'posterior_median', 15% 'prior_mode' or 16% 'prior mean'. 17% [double] np*1 array, values of the estimated parameters. 18% + controlled_varexo [cell] m*1 cell of row char array, list of controlled exogenous variables. 19% + conf_sig [double] scalar in [0,1], probability mass covered by the confidence bands. 20% 21% OUTPUTS 22% None. 23% 24% SPECIAL REQUIREMENTS 25% This routine has to be called after an estimation statement or an estimated_params block. 26% 27% REMARKS 28% [1] Results are stored in oo_.conditional_forecast. 29% [2] Use the function plot_icforecast to plot the results. 30 31% Copyright (C) 2006-2020 Dynare Team 32% 33% This file is part of Dynare. 34% 35% Dynare is free software: you can redistribute it and/or modify 36% it under the terms of the GNU General Public License as published by 37% the Free Software Foundation, either version 3 of the License, or 38% (at your option) any later version. 39% 40% Dynare is distributed in the hope that it will be useful, 41% but WITHOUT ANY WARRANTY; without even the implied warranty of 42% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 43% GNU General Public License for more details. 44% 45% You should have received a copy of the GNU General Public License 46% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 47 48global options_ oo_ M_ bayestopt_ estim_params_ 49 50if ~isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set) 51 if isfield(oo_,'posterior_mode') 52 options_cond_fcst.parameter_set = 'posterior_mode'; 53 elseif isfield(oo_,'mle_mode') 54 options_cond_fcst.parameter_set = 'mle_mode'; 55 else 56 error('No valid parameter set found') 57 end 58end 59 60if ~isfield(options_cond_fcst,'replic') || isempty(options_cond_fcst.replic) 61 options_cond_fcst.replic = 5000; 62end 63 64if ~isfield(options_cond_fcst,'periods') || isempty(options_cond_fcst.periods) 65 options_cond_fcst.periods = 40; 66end 67 68if ~isfield(options_cond_fcst,'conditional_forecast') || ~isfield(options_cond_fcst.conditional_forecast,'conf_sig') || isempty(options_cond_fcst.conditional_forecast.conf_sig) 69 options_cond_fcst.conditional_forecast.conf_sig = .8; 70end 71 72if isequal(options_cond_fcst.parameter_set,'calibration') 73 estimated_model = 0; 74else 75 estimated_model = 1; 76end 77 78if estimated_model 79 if options_.prefilter 80 error('imcforecast:: Conditional forecasting does not support the prefiltering option') 81 end 82 if ischar(options_cond_fcst.parameter_set) 83 switch options_cond_fcst.parameter_set 84 case 'posterior_mode' 85 xparam = get_posterior_parameters('mode',M_,estim_params_,oo_,options_); 86 graph_title='Posterior Mode'; 87 case 'posterior_mean' 88 xparam = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); 89 graph_title='Posterior Mean'; 90 case 'posterior_median' 91 xparam = get_posterior_parameters('median',M_,estim_params_,oo_,options_); 92 graph_title='Posterior Median'; 93 case 'mle_mode' 94 xparam = get_posterior_parameters('mode',M_,estim_params_,oo_,options_,'mle_'); 95 graph_title='ML Mode'; 96 case 'prior_mode' 97 xparam = bayestopt_.p5(:); 98 graph_title='Prior Mode'; 99 case 'prior_mean' 100 xparam = bayestopt_.p1; 101 graph_title='Prior Mean'; 102 otherwise 103 disp('imcforecast:: If the input argument is a string, then it has to be equal to:') 104 disp(' ''calibration'', ') 105 disp(' ''posterior_mode'', ') 106 disp(' ''posterior_mean'', ') 107 disp(' ''posterior_median'', ') 108 disp(' ''prior_mode'' or') 109 disp(' ''prior_mean''.') 110 error('imcforecast:: Wrong argument type!') 111 end 112 else 113 xparam = options_cond_fcst.parameter_set; 114 if length(xparam)~=length(M_.params) 115 error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!') 116 end 117 end 118 set_parameters(xparam); 119 [dataset_,dataset_info] = makedataset(options_); 120 data = transpose(dataset_.data); 121 data_index = dataset_info.missing.aindex; 122 gend = dataset_.nobs; 123 missing_value = dataset_info.missing.state; 124 125 %store qz_criterium 126 qz_criterium_old=options_.qz_criterium; 127 options_=select_qz_criterium_value(options_); 128 options_smoothed_state_uncertainty_old = options_.smoothed_state_uncertainty; 129 [atT, ~, ~, ~,ys, ~, ~, ~, ~, ~, ~, ~, ~, ~,M_,oo_,options_,bayestopt_] = ... 130 DsgeSmoother(xparam, gend, data, data_index, missing_value, M_, oo_, options_, bayestopt_, estim_params_); 131 options_.smoothed_state_uncertainty = options_smoothed_state_uncertainty_old; 132 %get constant part 133 if options_.noconstant 134 constant = zeros(size(ys,1),options_cond_fcst.periods+1); 135 else 136 if options_.loglinear 137 constant = repmat(log(ys),1,options_cond_fcst.periods+1); 138 else 139 constant = repmat(ys,1,options_cond_fcst.periods+1); 140 end 141 end 142 %get trend part (which also takes care of prefiltering); needs to 143 %include the last period 144 if bayestopt_.with_trend == 1 145 [trend_addition] =compute_trend_coefficients(M_,options_,size(bayestopt_.smoother_mf,1),gend+options_cond_fcst.periods); 146 trend_addition = trend_addition(:,gend:end); 147 else 148 trend_addition=zeros(size(bayestopt_.smoother_mf,1),1+options_cond_fcst.periods); 149 end 150 % add trend to constant 151 for obs_iter=1:length(options_.varobs) 152 j = strcmp(options_.varobs{obs_iter}, M_.endo_names); 153 constant(j,:) = constant(j,:) + trend_addition(obs_iter,:); 154 end 155 trend = constant(oo_.dr.order_var,:); 156 InitState(:,1) = atT(:,end); 157else 158 qz_criterium_old=options_.qz_criterium; 159 if isempty(options_.qz_criterium) 160 options_.qz_criterium = 1+1e-6; 161 end 162 graph_title='Calibration'; 163 if ~isfield(oo_.dr,'kstate') 164 error('You need to call stoch_simul before conditional_forecast') 165 end 166end 167 168if options_.logged_steady_state %if steady state was previously logged, undo this 169 oo_.dr.ys=exp(oo_.dr.ys); 170 oo_.steady_state=exp(oo_.steady_state); 171 options_.logged_steady_state=0; 172end 173 174[T, R, ys, ~, M_, options_, oo_] = dynare_resolve(M_, options_, oo_); 175 176if options_.loglinear && isfield(oo_.dr,'ys') && options_.logged_steady_state==0 %log steady state 177 oo_.dr.ys=log_variable(1:M_.endo_nbr,oo_.dr.ys,M_); 178 ys=oo_.dr.ys; 179 oo_.steady_state=log_variable(1:M_.endo_nbr,oo_.steady_state,M_); 180 options_.logged_steady_state=1; %set option for use in stoch_simul 181end 182 183if ~isdiag(M_.Sigma_e) 184 warning(sprintf('The innovations are correlated (the covariance matrix has non zero off diagonal elements), the results of the conditional forecasts will\ndepend on the ordering of the innovations (as declared after varexo) because a Cholesky decomposition is used to factorize the covariance matrix.\n\n=> It is preferable to declare the correlations in the model block (explicitly imposing the identification restrictions), unless you are satisfied\nwith the implicit identification restrictions implied by the Cholesky decomposition.')) 185 sQ = chol(M_.Sigma_e,'lower'); 186else 187 sQ = sqrt(M_.Sigma_e); 188end 189 190if ~estimated_model 191 if isempty(M_.endo_histval) 192 y0 = ys; 193 else 194 if options_.loglinear 195 %make sure that only states are updated (controls have value of 0 in vector) 196 y0=zeros(size(ys)); 197 y0_logged = log_variable(1:M_.endo_nbr,M_.endo_histval,M_); 198 y0(M_.endo_histval~=0)=y0_logged(M_.endo_histval~=0); 199 else 200 y0 = M_.endo_histval; 201 end 202 end 203 InitState(:,1) = y0(oo_.dr.order_var)-ys(oo_.dr.order_var,:); %initial state in deviations from steady state 204 trend = repmat(ys(oo_.dr.order_var,:),1,options_cond_fcst.periods+1); %trend needs to contain correct steady state 205end 206 207NumberOfStates = length(InitState); 208FORCS1 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic); 209 210FORCS1(:,1,:) = repmat(InitState,1,options_cond_fcst.replic); %set initial steady state to deviations from steady state in first period 211 212EndoSize = M_.endo_nbr; 213ExoSize = M_.exo_nbr; 214 215n1 = size(constrained_vars,1); 216n2 = size(options_cond_fcst.controlled_varexo,1); 217 218constrained_vars = oo_.dr.inv_order_var(constrained_vars); % must be in decision rule order 219 220if n1 ~= n2 221 error('imcforecast:: The number of constrained variables doesn''t match the number of controlled shocks') 222end 223 224% Get indices of controlled varexo. 225[~, controlled_varexo] = ismember(options_cond_fcst.controlled_varexo,M_.exo_names); 226 227mv = zeros(n1, NumberOfStates); 228mu = zeros(ExoSize, n2); 229for i=1:n1 230 mv(i,constrained_vars(i)) = 1; 231 mu(controlled_varexo(i),i) = 1; 232end 233 234% number of periods with constrained values 235cL = size(constrained_paths,2); 236 237%transform constrained periods into deviations from steady state; note that 238%trend includes last actual data point and therefore we need to start in 239%period 2 240constrained_paths = bsxfun(@minus,constrained_paths,trend(constrained_vars,2:1+cL)); 241 242FORCS1_shocks = zeros(n1,cL,options_cond_fcst.replic); 243 244%randn('state',0); 245 246for b=1:options_cond_fcst.replic %conditional forecast using cL set to constrained values 247 shocks = sQ*randn(ExoSize,options_cond_fcst.periods); 248 shocks(controlled_varexo,:) = zeros(n1, options_cond_fcst.periods); 249 [FORCS1(:,:,b), FORCS1_shocks(:,:,b)] = mcforecast3(cL,options_cond_fcst.periods,constrained_paths,shocks,FORCS1(:,:,b),T,R,mv, mu); 250 FORCS1(:,:,b)=FORCS1(:,:,b)+trend; %add trend 251end 252if max(max(max(abs(bsxfun(@minus,FORCS1(constrained_vars,1+1:1+cL,:),trend(constrained_vars,1:cL)+constrained_paths)))))>1e-4 253 fprintf('\nconditional_forecasts: controlling of variables was not successful.\n') 254 fprintf('This can be due to numerical imprecision (e.g. explosive simulations)\n') 255 fprintf('or because the instrument(s) do not allow controlling the variable(s).\n') 256end 257mFORCS1 = mean(FORCS1,3); 258mFORCS1_shocks = mean(FORCS1_shocks,3); 259 260tt = (1-options_cond_fcst.conditional_forecast.conf_sig)/2; 261t1 = max(1,round(options_cond_fcst.replic*tt)); 262t2 = min(options_cond_fcst.replic,round(options_cond_fcst.replic*(1-tt))); 263 264forecasts.controlled_variables = constrained_vars; 265forecasts.instruments = options_cond_fcst.controlled_varexo; 266 267for i = 1:EndoSize 268 forecasts.cond.Mean.(M_.endo_names{oo_.dr.order_var(i)}) = mFORCS1(i,:)'; 269 if size(FORCS1,2)>1 270 tmp = sort(squeeze(FORCS1(i,:,:))'); 271 else 272 tmp = sort(squeeze(FORCS1(i,:,:))); 273 end 274 forecasts.cond.ci.(M_.endo_names{oo_.dr.order_var(i)}) = [tmp(t1,:)' ,tmp(t2,:)' ]'; 275end 276 277for i = 1:n1 278 forecasts.controlled_exo_variables.Mean.(options_cond_fcst.controlled_varexo{i}) = mFORCS1_shocks(i,:)'; 279 if size(FORCS1_shocks,2)>1 280 tmp = sort(squeeze(FORCS1_shocks(i,:,:))'); 281 else 282 tmp = sort(squeeze(FORCS1_shocks(i,:,:))); 283 end 284 forecasts.controlled_exo_variables.ci.(options_cond_fcst.controlled_varexo{i}) = [tmp(t1,:)' ,tmp(t2,:)' ]'; 285end 286 287clear FORCS1 mFORCS1_shocks; 288 289FORCS2 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic); 290FORCS2(:,1,:) = repmat(InitState,1,options_cond_fcst.replic); %set initial steady state to deviations from steady state in first period 291 292for b=1:options_cond_fcst.replic %conditional forecast using cL set to 0 293 shocks = sQ*randn(ExoSize,options_cond_fcst.periods); 294 shocks(controlled_varexo,:) = zeros(n1, options_cond_fcst.periods); 295 FORCS2(:,:,b) = mcforecast3(0,options_cond_fcst.periods,constrained_paths,shocks,FORCS2(:,:,b),T,R,mv, mu)+trend; 296end 297 298mFORCS2 = mean(FORCS2,3); 299 300for i = 1:EndoSize 301 forecasts.uncond.Mean.(M_.endo_names{oo_.dr.order_var(i)})= mFORCS2(i,:)'; 302 if size(FORCS2,2)>1 303 tmp = sort(squeeze(FORCS2(i,:,:))'); 304 else 305 tmp = sort(squeeze(FORCS2(i,:,:))); 306 end 307 forecasts.uncond.ci.(M_.endo_names{oo_.dr.order_var(i)}) = [tmp(t1,:)' ,tmp(t2,:)' ]'; 308end 309forecasts.graph.title = graph_title; 310forecasts.graph.fname = M_.fname; 311 312%reset qz_criterium 313options_.qz_criterium=qz_criterium_old; 314oo_.conditional_forecast = forecasts; 315