diff --git a/doc/dynare.texi b/doc/dynare.texi
index 5381c24a7feb366dfcd4e0c7135f01bb2721ad08..041cebd7e5e45138e323e23df39f3cc8084d2865 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -2404,10 +2404,6 @@ values (xx);
end;
@end example
-In case of conditional forecasts using an extended path method,
-the shock command is extended in order to determine the nature of the expectation
-on the endogenous variable paths. The command @code{expectation} has, in this case,
-to be added after the @code{values} command (@pxref{Forecasting}).
@customhead{In stochastic context}
@@ -4210,7 +4206,7 @@ apply.
The following line:
@example
-corr eps_1, eps_2, 0.5, , , beta_pdf, 0, 0.3, -1, 1;
+corr eps_1, eps_2, 0.5, , , beta_pdf, 0, 0.3, -1, 1;
@end example
sets a generalized beta prior for the correlation between @code{eps_1} and
@code{eps_2} with mean 0 and variance 0.3. By setting
@@ -4222,7 +4218,7 @@ respectively. The initial value is set to 0.5.
Similarly, the following line:
@example
-corr eps_1, eps_2, 0.5, -0.5, 1, beta_pdf, 0, 0.3, -1, 1;
+corr eps_1, eps_2, 0.5, -0.5, 1, beta_pdf, 0, 0.3, -1, 1;
@end example
sets the same generalized beta distribution as before, but now truncates
this distribution to [-0.5,1] through the use of @var{LOWER_BOUND} and
@@ -5559,13 +5555,6 @@ DSGE model, by finding the structural shocks that are needed to match
the restricted paths. Use @code{conditional_forecast},
@code{conditional_forecast_paths} and @code{plot_conditional_forecast}
for that purpose.
-If the model contains strong non-linearities, the conditional forecasts
-can be computed using an extended path method with the @code{simulation_type}
-option in @code{conditional_forecast} command set to @code{deterministic}.
-Because in this case deterministic simulations are carried out,
-the nature of the shocks (surprise or perfect foresight) has to be indicated
-in the @code{conditional_forecast_paths} block, using the command @code{expectation}
-for each endogenous path. The forecasts are plotted using the @code{rplot} command.
Finally, it is possible to do forecasting with a Bayesian VAR using
the @code{bvar_forecast} command.
@@ -5726,10 +5715,7 @@ command has to be called after estimation.
Use @code{conditional_forecast_paths} block to give the list of
constrained endogenous, and their constrained future path.
-If an extended path method is applied on the original dsge model,
-the nature of the expectation on the constrained endogenous has to be
-specified using expectation command. Option
-@code{controlled_varexo} is used to specify the structural shocks
+Option @code{controlled_varexo} is used to specify the structural shocks
which will be matched to generate the constrained path.
Use @code{plot_conditional_forecast} to graph the results.
@@ -5756,14 +5742,6 @@ Number of simulations. Default: @code{5000}.
@item conf_sig = @var{DOUBLE}
Level of significance for confidence interval. Default: @code{0.80}
-@item simulation_type = @code{stochastic} | @code{deterministic}
-Indicates the nature of simulations used to compute the conditional forecast.
-The default value @code{stochastic} is used, when simulations are computed
-using the reduced form representation of the DSGE model.
-If the model has to be simulated using extended path method on the original
-DSGE model, @code{simulation_type} has to be set equal to @code{deterministic}.
-
-
@end table
@outputhead
@@ -5833,31 +5811,6 @@ plot_conditional_forecast(periods = 10) a y;
@end example
-@examplehead
-@example
-/* conditional forecast using extended path method
-with perfect foresight on r path*/
-var y r
-varexo e u;
-
-@dots{}
-
-conditional_forecast_paths;
-var y;
-periods 1:3, 4:5;
-values 2, 5;
-var r
-periods 1:5;
-values 3;
-expectation perfect_foresight;
-end;
-
-conditional_forecast(parameter_set = calibration, controlled_varexo = (e, u), simulation_type=deterministic);
-
-rplot a;
-rplot y;
-@end example
-
@end deffn
@deffn Block conditional_forecast_paths ;
@@ -5870,10 +5823,6 @@ example.
The syntax of the block is the same than the deterministic shocks in
the @code{shocks} blocks (@pxref{Shocks on exogenous variables}).
-If the conditional forecast is carried out using the extended path method
-on the original DSGE model, the nature of the expectation have to be specified
-for each endogenous path, using the @code{expectation} = @code{surprise} | @code{perfect_foresight}.
-By default, @code{expectation} is equal to @code{surprise}.
@end deffn
@@ -5907,6 +5856,67 @@ See @file{bvar-a-la-sims.pdf}, which comes with Dynare distribution,
for more information on this command.
@end deffn
+If the model contains strong non-linearities or if some perfectly expected shocks are considered, the forecasts and the conditional forecasts
+can be computed using an extended path method. The forecast scenario describing the shocks and/or the constrained paths on some endogenous variables should be build.
+The first step is the forecast scenario initialization using the function @code{init_plan}:
+
+@anchor{init_plan}
+@deffn {MATLAB/Octave command} HANDLE = init_plan(DATES) ;
+
+Creates a new forecast scenario for a forecast period (indicated as a dates class, see @ref{dates class members}). This function return a handle on the new forecast scenario.
+
+@end deffn
+
+The forecast scenario can contain some simple shocks on the exogenous variables. This shocks are described using the function @code{basic_plan}:
+
+@anchor{basic_plan}
+@deffn {MATLAB/Octave command} HANDLE = basic_plan(HANDLE, 'VARIABLE_NAME', 'SHOCK_TYPE', DATES, MATLAB VECTOR OF DOUBLE | [DOUBLE | EXPRESSION [DOUBLE | | EXPRESSION] ] ) ;
+
+Adds to the forecast scenario a shock on the exogenous variable indicated between quotes in the second argument. The shock type has to be specified in the third argument between quotes: 'surprise' in case of an unexpected shock or 'perfect_foresight' for a perfectly anticipated shock. The fourth argument indicates the period of the shock using a dates class (see @ref{dates class members}). The last argument is the shock path indicated as a Matlab vector of double. This function return the handle of the updated forecast scenario.
+
+@end deffn
+
+The forecast scenario can also contain a constrained path on an endogenous variable. The values of the related exogenous variable compatible with the constrained path are in this case computed. In other words, a conditional forecast is performed. This kind of shock is described with the function @code{flip_plan}:
+
+@anchor{flip_plan}
+@deffn {MATLAB/Octave command} HANDLE = flip_plan(HANDLE, 'VARIABLE_NAME, 'VARIABLE_NAME', 'SHOCK_TYPE', DATES, MATLAB VECTOR OF DOUBLE | [DOUBLE | EXPRESSION [DOUBLE | | EXPRESSION] ] ) ;
+
+Adds to the forecast scenario a constrained path on the endogenous variable specified between quotes in the second argument. The associated exogenous variable provided in the third argument between quotes, is considered as an endogenous variable and its values compatible with the constrained path on the endogenous variable will be computed. The nature of the expectation on the constrained path has to be specified in the fourth argument between quotes: 'surprise' in case of an unexpected path or 'perfect_foresight' for a perfectly anticipated path. The fifth argument indicates the period where the path of the endogenous variable is constrained using a dates class (see @ref{dates class members}). The last argument contains the constrained path as a Matlab vector of double. This function return the handle of the updated forecast scenario.
+
+@end deffn
+
+Once the forecast scenario if fully described, the forecast is computed with the command @code{det_cond_forecast}:
+@anchor{det_cond_forecast}
+@deffn {MATLAB/Octave command} DSERIES = det_cond_forecast(HANDLE[, DSERIES [, DATES]]) ;
+
+Computes the forecast or the conditional forecast using an extended path method for the given forecast scenario (first argument). The past values of the endogenous and exogenous variables provided with a dseries class (see @ref{dseries class members}) can be indicated in the second argument. By default, the past values of the variables are equal to their steady-state values. The initial date of the forecast can be provided in the third argument. By default, the forecast will start at the first date indicated in the @code{init_plan} command. This function returns a dset containing the historical and forecast values for the endogenous and exogenous variables.
+
+@end deffn
+
+
+
+@examplehead
+@example
+/* conditional forecast using extended path method
+with perfect foresight on r path*/
+var y r
+varexo e u;
+
+@dots{}
+
+smoothed = dseries('smoothed_variables.csv');
+
+fplan = init_plan(2013Q4:2029Q4);
+
+fplan = flip_plan(fplan, 'y', 'u', 'surprise', 2013Q4:2014Q4, [1 1.1 1.2 1.1 ]);
+
+fplan = flip_plan(fplan, 'r', 'e', 'perfect_foresight', 2013Q4:2014Q4, [2 1.9 1.9 1.9 ]);
+
+dset_forecast = det_cond_forecast(fplan, smoothed);
+
+plot(dset_forecast.{'y','u'});
+plot(dset_forecast.{'r','e'});
+@end example
@node Optimal policy
@section Optimal policy
diff --git a/matlab/basic_plan.m b/matlab/basic_plan.m
new file mode 100644
index 0000000000000000000000000000000000000000..a2e11477f4cbd48e34012a24ca0e498e6527a3a3
--- /dev/null
+++ b/matlab/basic_plan.m
@@ -0,0 +1,82 @@
+function plan = basic_plan(plan, exogenous, expectation_type, date, value)
+% Adds a simple shock to the forecast scenario plan
+%
+% INPUTS
+% o plan [structure] A structure describing the different shocks and the implied variables, the date of the shocks and the path of the shock (forecast scenario).
+% The plan structure is created by the functions init_plan, basic_plan and flip_plan
+% o exogenous [string] A string containing the name of the exognous shock.
+% o expectation_type [string] A string indicating the type of expectation: 'surprise' for an unexpected shock, and 'perfect_foresight' for a perfectly anticpated shock.
+% o date [dates] The period of the shock
+% o value [array of double] A vector of double containing the values of the exogenous variable
+%
+%
+% OUTPUTS
+% plan [structure] Returns a structure containing the updated forecast scenario.
+%
+%
+% Copyright (C) 2013-2014 Dynare Team
+%
+% 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/>.
+
+ if ~ischar(expectation_type) || size(expectation_type,1) ~= 1
+ error(['in basic_plan the third argument should be a string containing the simulation type (''perfect_foresight'' or ''surprise'')']);
+ end
+ exogenous = strtrim(exogenous);
+ ix = find(strcmp(exogenous, plan.exo_names));
+ if isempty(ix)
+ error(['in basic_plan the second argument ' exogenous ' is not an exogenous variable']);
+ end;
+ sdate = length(date);
+ if sdate > 1
+ if date(1) < plan.date(1) || date(end) > plan.date(end)
+ error(['in basic_plan the forth argument (date=' date ') must lay inside the plan.date ' plan.date]);
+ end
+ else
+ if date < plan.date(1) || date > plan.date(end)
+ error(['in basic_plan the forth argument (date=' date ') must lay iside the plan.date ' plan.date]);
+ end
+ end
+ if ~isempty(plan.options_cond_fcst_.controlled_varexo)
+ common_var = find(ix == plan.options_cond_fcst_.controlled_varexo);
+ if ~isempty(common_var)
+ common_date = intersect(date, plan.constrained_date_{common_var});
+ if ~isempty(common_date)
+ if common_date.length > 1
+ the_dates = [cell2mat(strings(common_date(1))) ':' cell2mat(strings(common_date(end)))];
+ else
+ the_dates = cell2mat(strings(common_date));
+ end
+ error(['Impossible case: ' plan.exo_names{plan.options_cond_fcst_.controlled_varexo(common_var)} ' is used both as a shock and as an endogenous variable to control the path of ' plan.endo_names{plan.constrained_vars_(common_var)} ' at the dates ' the_dates]);
+ end
+ end
+ end
+ if isempty(plan.shock_vars_)
+ plan.shock_vars_ = ix;
+ if strcmp(expectation_type, 'perfect_foresight')
+ plan.shock_perfect_foresight_ = 1;
+ else
+ plan.shock_perfect_foresight_ = 0;
+ end
+ else
+ plan.shock_vars_ = [plan.shock_vars_ ; ix];
+ if strcmp(expectation_type, 'perfect_foresight')
+ plan.shock_perfect_foresight_ = [plan.shock_perfect_foresight_ ; 1];
+ else
+ plan.shock_perfect_foresight_ = [plan.shock_perfect_foresight_ ; 0];
+ end
+ end
+ plan.shock_date_{length(plan.shock_date_) + 1} = date;
+ plan.shock_paths_{length(plan.shock_paths_) + 1} = value;
diff --git a/matlab/det_cond_forecast.m b/matlab/det_cond_forecast.m
index 9cfc77cbba0f3c776da35b4164be28e774761d12..4cbbe125fe764ec249899559553b3d0015e21d29 100644
--- a/matlab/det_cond_forecast.m
+++ b/matlab/det_cond_forecast.m
@@ -1,357 +1,846 @@
-function det_cond_forecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight)
-% Computes conditional forecasts for a deterministic model.
-%
-% INPUTS
-% o constrained_paths [double] m*p array, where m is the number of constrained endogenous variables and p is the number of constrained periods.
-% o constrained_vars [char] m*x array holding the names of the controlled endogenous variables.
-% o options_cond_fcst [structure] containing the options. The fields are:
-% + replic [integer] scalar, number of monte carlo simulations.
-% + parameter_set [char] values of the estimated parameters:
-% "posterior_mode",
-% "posterior_mean",
-% "posterior_median",
-% "prior_mode" or
-% "prior mean".
-% [double] np*1 array, values of the estimated parameters.
-% + controlled_varexo [char] m*x array, list of controlled exogenous variables.
-% + conf_sig [double] scalar in [0,1], probability mass covered by the confidence bands.
-% o constrained_perfect_foresight [double] m*1 array indicating if the endogenous variables path is perfectly foresight (1) or is a surprise (0)
-%
-%
-% OUTPUTS
-% None.
-%
-% SPECIAL REQUIREMENTS
-% This routine has to be called after an estimation statement or an estimated_params block.
-%
-% REMARKS
-% [1] Results are stored in a structure which is saved in a mat file called conditional_forecasts.mat.
-% [2] Use the function plot_icforecast to plot the results.
-
-% Copyright (C) 2013 Dynare Team
-%
-% 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/>.
-
-global options_ oo_ M_
-
-if ~isfield(options_cond_fcst,'periods') || isempty(options_cond_fcst.periods)
- options_cond_fcst.periods = 100;
-end
-
-maximum_lag = M_.maximum_lag;
-maximum_lead = M_.maximum_lead;
-ys = oo_.steady_state;
-ny = size(ys,1);
-xs = [oo_.exo_steady_state ; oo_.exo_det_steady_state];
-nx = size(xs,1);
-
-constrained_periods = size(constrained_paths,2);
-n_endo_constrained = size(constrained_vars,1);
-if isfield(options_cond_fcst,'controlled_varexo')
- n_control_exo = size(options_cond_fcst.controlled_varexo, 1);
- if n_control_exo ~= n_endo_constrained
- error(['det_cond_forecast:: the number of exogenous controlled variables (' int2str(n_control_exo) ') has to be equal to the number of constrained endogenous variabes (' int2str(n_endo_constrained) ')'])
- end;
-else
- error('det_cond_forecast:: to run a deterministic conditional forecast you have to specified the exogenous variables controlled using the option controlled_varex in forecast command');
-end;
-
-exo_names = M_.exo_names;
-controlled_varexo = zeros(1,n_control_exo);
-for i = 1:nx
- for j=1:n_control_exo
- if strcmp(deblank(exo_names(i,:)), deblank(options_cond_fcst.controlled_varexo(j,:)))
- controlled_varexo(j) = i;
- end
- end
-end
-
-%todo check if zero => error message
-
-save_options_initval_file = options_.initval_file;
-options_.initval_file = '__';
-
-[pos_constrained_pf, junk] = find(constrained_perfect_foresight);
-indx_endo_solve_pf = constrained_vars(pos_constrained_pf);
-if isempty(indx_endo_solve_pf)
- pf = 0;
-else
- pf = length(indx_endo_solve_pf);
-end;
-indx_endo_solve_surprise = setdiff(constrained_vars, indx_endo_solve_pf);
-
-if isempty(indx_endo_solve_surprise)
- surprise = 0;
-else
- surprise = length(indx_endo_solve_surprise);
-end;
-
-eps = options_.solve_tolf;
-maxit = options_.simul.maxit;
-
-initial_conditions = oo_.steady_state;
-past_val = 0;
-save_options_periods = options_.periods;
-options_.periods = options_cond_fcst.periods;
-save_options_dynatol_f = options_.dynatol.f;
-options_.dynatol.f = 1e-8;
-eps1 = 1e-4;
-exo = zeros(maximum_lag + options_cond_fcst.periods, nx);
-endo = zeros(maximum_lag + options_cond_fcst.periods, ny);
-endo(1,:) = oo_.steady_state';
-% if all the endogenous paths are perfectly anticipated we do not need to
-% implement extended path
-if pf && ~surprise
- time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods;
- second_system_size = pf * constrained_periods;
- J = zeros(second_system_size,second_system_size);
- r = zeros(second_system_size,1);
- indx_endo = zeros(second_system_size,1);
- col_count = 1;
- for j = 1:length(constrained_vars)
- indx_endo(col_count : col_count + constrained_periods - 1) = constrained_vars(j) + (time_index_constraint - 1) * ny;
- col_count = col_count + constrained_periods;
- end;
- make_ex_;
- make_y_;
- it = 1;
- convg = 0;
- while ~convg && it <= maxit
- disp('---------------------------------------------------------------------------------------------');
- disp(['iteration ' int2str(it)]);
- not_achieved = 1;
- alpha = 1;
- while not_achieved
- simul();
- result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * constrained_periods;
- if (~oo_.deterministic_simulation.status || ~result) && it > 1
- not_achieved = 1;
- alpha = alpha / 2;
- col_count = 1;
- for j = controlled_varexo
- oo_.exo_simul(time_index_constraint,j) = (old_exo(:,j) + alpha * D_exo(col_count: col_count + constrained_periods - 1));
- col_count = col_count + constrained_periods;
- end;
- disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha)]);
- oo_.endo_simul = repmat(oo_.steady_state, 1, options_cond_fcst.periods + 2);
- else
- not_achieved = 0;
- end;
- end;
-
- y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
- ys = oo_.endo_simul(indx_endo);
- col_count = 1;
- for j = 1:length(constrained_vars)
- y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
- r(col_count:col_count+constrained_periods - 1) = (y - constrained_paths(j,1:constrained_periods))';
- col_count = col_count + constrained_periods;
- end;
- col_count = 1;
- for j = controlled_varexo
- for time = time_index_constraint
- saved = oo_.exo_simul(time,j);
- oo_.exo_simul(time,j) = oo_.exo_simul(time,j) + eps1;
- simul();
- J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
- oo_.exo_simul(time,j) = saved;
- col_count = col_count + 1;
- end;
- end;
- normr = norm(r, 1);
-
- disp(['iteration ' int2str(it) ' error = ' num2str(normr)]);
-
- if normr <= eps
- convg = 1;
- disp('convergence achieved');
- else
- % Newton update on exogenous shocks
- old_exo = oo_.exo_simul(time_index_constraint,:);
- D_exo = - J \ r;
- col_count = 1;
- for j = controlled_varexo
- oo_.exo_simul(time_index_constraint,j) = (oo_.exo_simul(time_index_constraint,j) + D_exo(col_count: col_count + constrained_periods - 1));
- col_count = col_count + constrained_periods;
- end;
- end;
- it = it + 1;
- end;
- endo(maximum_lag + 1 :maximum_lag + options_cond_fcst.periods ,:) = oo_.endo_simul(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods )';
- exo = oo_.exo_simul;
-else
- for t = 1:constrained_periods
- disp('=============================================================================================');
- disp(['t=' int2str(t) ' constrained_paths(:,t)=' num2str(constrained_paths(:,t)') ]);
- disp('=============================================================================================');
- if t == 1
- make_ex_;
- exo_init = oo_.exo_simul;
- make_y_;
- end;
- oo_.exo_simul = exo_init;
- oo_.endo_simul(:,1) = initial_conditions;
- convg = 0;
- it = 1;
-
- time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods - t + 1;
-
- second_system_size = surprise + pf * (constrained_periods - t + 1);
- indx_endo = zeros(second_system_size,1);
- col_count = 1;
- for j = 1:length(constrained_vars)
- if constrained_perfect_foresight(j)
- indx_endo(col_count : col_count + constrained_periods - t) = constrained_vars(j) + (time_index_constraint - 1) * ny;
- col_count = col_count + constrained_periods - t + 1;
- else
- indx_endo(col_count) = constrained_vars(j) + maximum_lag * ny;
- col_count = col_count + 1;
- end;
- end;
-
- J = zeros(second_system_size,second_system_size);
-
- r = zeros(second_system_size,1);
-
- while ~convg && it <= maxit
- disp('---------------------------------------------------------------------------------------------');
- disp(['iteration ' int2str(it) ' time ' int2str(t)]);
- not_achieved = 1;
- alpha = 1;
- while not_achieved
- simul();
- result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * length(time_index_constraint);
- if (~oo_.deterministic_simulation.status || ~result) && it > 1
- not_achieved = 1;
- alpha = alpha / 2;
- col_count = 1;
- for j = controlled_varexo
- if constrained_perfect_foresight(j)
- oo_.exo_simul(time_index_constraint,j) = (old_exo(time_index_constraint,j) + alpha * D_exo(col_count: col_count + constrained_periods - t));
- col_count = col_count + constrained_periods - t + 1;
- else
- oo_.exo_simul(maximum_lag + 1,j) = old_exo(maximum_lag + 1,j) + alpha * D_exo(col_count);
- col_count = col_count + 1;
- end;
- end;
- disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha) ' result=' num2str(result) ' sum(sum)=' num2str(sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint))))) ' ny * length(time_index_constraint)=' num2str(ny * length(time_index_constraint)) ' oo_.deterministic_simulation.status=' num2str(oo_.deterministic_simulation.status)]);
- oo_.endo_simul = [initial_conditions repmat(oo_.steady_state, 1, options_cond_fcst.periods + 1)];
- else
- not_achieved = 0;
- end;
- end;
- if t==constrained_periods
- ycc = oo_.endo_simul;
- end;
- yc = oo_.endo_simul(:,maximum_lag + 1);
- ys = oo_.endo_simul(indx_endo);
-
- col_count = 1;
- for j = 1:length(constrained_vars)
- if constrained_perfect_foresight(j)
- y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
- r(col_count:col_count+constrained_periods - t) = (y - constrained_paths(j,t:constrained_periods))';
- col_count = col_count + constrained_periods - t + 1;
- else
- y = yc(constrained_vars(j));
- r(col_count) = y - constrained_paths(j,t);
- col_count = col_count + 1;
- end;
- end
-
- disp('computation of derivatives w.r. to exogenous shocks');
- col_count = 1;
- for j = 1:length(controlled_varexo)
- j_pos = controlled_varexo(j);
- if constrained_perfect_foresight(j)
- for time = time_index_constraint
- saved = oo_.exo_simul(time,j_pos);
- oo_.exo_simul(time,j_pos) = oo_.exo_simul(time,j_pos) + eps1;
- simul();
- J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
- oo_.exo_simul(time,j_pos) = saved;
- col_count = col_count + 1;
- end;
- else
- saved = oo_.exo_simul(maximum_lag+1,j_pos);
- oo_.exo_simul(maximum_lag+1,j_pos) = oo_.exo_simul(maximum_lag+1,j_pos) + eps1;
- simul();
- J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
- oo_.exo_simul(maximum_lag+1,j_pos) = saved;
- col_count = col_count + 1;
- end;
- end;
-
- normr = norm(r, 1);
-
- disp(['iteration ' int2str(it) ' error = ' num2str(normr) ' at time ' int2str(t)]);
-
- if normr <= eps
- convg = 1;
- disp('convergence achieved');
- else
- % Newton update on exogenous shocks
- D_exo = - J \ r;
- old_exo = oo_.exo_simul;
- col_count = 1;
- for j = 1:length(controlled_varexo)
- j_pos=controlled_varexo(j);
- if constrained_perfect_foresight(j)
- oo_.exo_simul(time_index_constraint,j_pos) = (oo_.exo_simul(time_index_constraint,j_pos) + D_exo(col_count: col_count + constrained_periods - t));
- col_count = col_count + constrained_periods - t + 1;
- else
- oo_.exo_simul(maximum_lag + 1,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos) + D_exo(col_count);
- col_count = col_count + 1;
- end;
- end;
- end;
- it = it + 1;
- end;
- if ~convg
- error(['convergence not achived at time ' int2str(t) ' after ' int2str(it) ' iterations']);
- end;
- for j = 1:length(controlled_varexo)
- j_pos=controlled_varexo(j);
- if constrained_perfect_foresight(j)
- % in case of mixed surprise and perfect foresight
- % endogenous path at each date all the exogenous paths have to be
- % stored. The paths are stacked in exo.
- for time = time_index_constraint;
- exo(past_val + time,j_pos) = oo_.exo_simul(time,j_pos);
- end
- else
- exo(maximum_lag + t,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos);
- end;
- end;
- past_val = past_val + length(time_index_constraint);
- if t < constrained_periods
- endo(maximum_lag + t,:) = yc;
- else
- endo(maximum_lag + t :maximum_lag + options_cond_fcst.periods ,:) = ycc(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods - constrained_periods + 1)';
- end;
- initial_conditions = yc;
- end;
-end;
-options_.periods = save_options_periods;
-options_.dynatol.f = save_options_dynatol_f;
-options_.initval_file = save_options_initval_file;
-% disp('endo');
-% disp(endo);
-% disp('exo');
-% disp(exo);
-
-oo_.endo_simul = endo';
-oo_.exo_simul = exo;
+function data_set = det_cond_forecast(varargin)
+% Computes conditional forecasts using the extended path method.
+%
+% INPUTS
+% o plan [structure] A structure describing the different shocks and the endogenous varibales, the date of the shocks and the path of the shock.
+% The plan structure is created by the functions init_plan, basic_plan and flip_plan
+% o [dataset] [dseries] A dserie containing the initial values of the shocks and the endogenous variables (usually the dseries generated by the smoother).
+% o [starting_date] [dates] The first date of the forecast.
+%
+%
+% OUTPUTS
+% dataset [dseries] Returns a dseries containing the forecasted endgenous variables and shocks
+%
+%
+% Copyright (C) 2013-2014 Dynare Team
+%
+% 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/>.
+
+global options_ oo_ M_
+pp = 2;
+initial_conditions = oo_.steady_state;
+
+%We have to get an initial guess for the conditional forecast
+% and terminal conditions for the non-stationary variables, we
+% use the first order approximation of the rational expectation solution.
+dr = struct();
+oo_.dr=set_state_space(dr,M_,options_);
+options_.order = 1;
+[dr,Info,M_,options_,oo_] = resol(0,M_,options_,oo_);
+b_surprise = 0;
+b_pf = 0;
+surprise = 0;
+pf = 0;
+is_shock = [];
+is_constraint = [];
+if length(varargin) > 3
+ % regular way to call
+ constrained_paths = varargin{1};
+ max_periods_simulation = size(constrained_paths, 2);
+ constrained_vars = varargin{2};
+ options_cond_fcst = varargin{3};
+ constrained_perfect_foresight = varargin{4};
+ constraint_index = cell(max_periods_simulation,1);
+ nvars = length(constrained_vars);
+ for i = 1:max_periods_simulation
+ constraint_index{i} = 1:nvars;
+ end;
+ direct_mode = 0;
+ shocks_present = 0;
+ controlled_varexo = options_cond_fcst.controlled_varexo;
+ nvarexo = size(controlled_varexo, 1);
+ options_cond_fcst.controlled_varexo = zeros(nvarexo,1);
+ exo_names = cell(M_.exo_nbr,1);
+ for i = 1:M_.exo_nbr
+ exo_names{i} = deblank(M_.exo_names(i,:));
+ end
+ for i = 1:nvarexo
+ j = find(strcmp(controlled_varexo(i,:), exo_names));
+ if ~isempty(j)
+ options_cond_fcst.controlled_varexo(i) = j;
+ else
+ error(['Unknown exogenous variable ' controlled_varexo(i,:)]);
+ end
+ end
+
+else
+ % alternative way to call
+ plan = varargin{1};
+ if length(varargin) >= 2
+ dset = varargin{2};
+ if ~isa(dset,'dseries')
+ error('the second argmuent should be a dseries');
+ end
+ if length(varargin) >= 3
+ range = varargin{3};
+ if ~isa(range,'dates')
+ error('the third argmuent should be a dates');
+ end
+ %if (range(range.ndat) > dset.time(dset.nobs) )
+ if (range(range.ndat) > dset.dates(dset.nobs) )
+ s1 = strings(dset.dates(dset.nobs));
+ s2 = strings(range(range.ndat));
+ error(['the dseries ' inputname(2) ' finish at time ' s1{1} ' before the last period of forecast ' s2{1}]);
+ end
+
+ sym_dset = dset(dates(-range(1)):dates(range(range.ndat)));
+
+ oo_.exo_simul = repmat(oo_.exo_steady_state',max(range.ndat + 1, options_.periods),1);
+ oo_.endo_simul = repmat(oo_.steady_state, 1, max(range.ndat + 1, options_.periods));
+
+ for i = 1:sym_dset.vobs
+ iy = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.endo_names)));
+ if ~isempty(iy)
+ oo_.endo_simul(iy,1:sym_dset.nobs) = sym_dset.data(:, i);
+ initial_conditions(iy) = sym_dset.data(1, i);
+ else
+ ix = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.exo_names)));
+ if ~isempty(ix)
+ oo_.exo_simul(1, ix) = sym_dset.data(1, i)';
+ else
+ %warning(['The variable ' sym_dset.name{i} ' in the dataset ' inputname(2) ' is not a endogenous neither an exogenous variable!!']);
+ end
+ end
+ end
+ for i = 1:length(M_.aux_vars)
+ if M_.aux_vars(i).type == 1 %lag variable
+ iy = find(strcmp(deblank(M_.endo_names(M_.aux_vars(i).orig_index,:)), sym_dset.name));
+ if ~isempty(iy)
+ oo_.endo_simul(M_.aux_vars(i).endo_index, 1:sym_dset.nobs) = dset(dates(range(1) + (M_.aux_vars(i).orig_lead_lag - 1)):dates(range(range.ndat) + M_.aux_vars(i).orig_lead_lag)).data(:,iy);
+ initial_conditions(M_.aux_vars(i).endo_index) = dset(dates(range(1) + (M_.aux_vars(i).orig_lead_lag - 1))).data(:,iy);
+ else
+ warning(['The variable auxiliary ' M_.endo_names(M_.aux_vars(i).endo_index, :) ' associated to the variable ' M_.endo_names(M_.aux_vars(i).orig_index,:) ' do not appear in the dataset']);
+ end
+ else
+ oo_.endo_simul(M_.aux_vars(i).endo_index, 1:sym_dset.nobs) = repmat(oo_.steady_state(M_.aux_vars(i).endo_index), 1, range.ndat + 1);
+ end
+ end
+ else
+ error('impossible case');
+ end;
+
+ else
+ oo_.exo_simul = repmat(oo_.exo_steady_state',options_.periods+2,1);
+ oo_.endo_simul = repmat(oo_.steady_state, 1, options_.periods+2);
+ end
+
+ direct_mode = 1;
+ constrained_paths = plan.constrained_paths_;
+ constrained_vars = plan.constrained_vars_;
+ options_cond_fcst = plan.options_cond_fcst_;
+ constrained_perfect_foresight = plan.constrained_perfect_foresight_;
+ constrained_periods = plan.constrained_date_;
+ if ~isempty(plan.shock_paths_)
+ shock_paths = plan.shock_paths_;
+ shock_vars = plan.shock_vars_;
+ shock_perfect_foresight = plan.shock_perfect_foresight_;
+ shock_periods = plan.shock_date_;
+ shocks_present = 1;
+ else
+ shocks_present = 0;
+ end
+
+ total_periods = plan.date;
+
+end;
+
+if ~isfield(options_cond_fcst,'periods') || isempty(options_cond_fcst.periods)
+ options_cond_fcst.periods = 100;
+end
+
+options_.periods = 10;
+
+if direct_mode == 1
+ n_periods = length(constrained_periods);
+ is_constraint = zeros(size(total_periods), n_periods);
+ constrained_paths_cell = constrained_paths;
+ clear constrained_paths;
+ constrained_paths = zeros(n_periods, size(total_periods));
+ max_periods_simulation = 0;
+ for i = 1:n_periods
+ period_i = constrained_periods{i};
+ %period_i
+ tp = total_periods(1);
+ if size(period_i) > 1
+ init_periods = period_i(1);
+ tp_end = period_i(end);
+ else
+ init_periods = period_i;
+ tp_end = period_i;
+ end;
+ tp0 = tp;
+ while tp < init_periods
+ tp = tp + 1;
+ end
+ j = 0;
+ while tp <= tp_end
+ is_constraint(tp - tp0 + 1, i) = 1;
+ constrained_paths(i, tp - tp0 + 1) = constrained_paths_cell{i}(j + 1);
+ tp = tp + 1;
+ j = j + 1;
+ end;
+ if tp - tp0 > max_periods_simulation
+ max_periods_simulation = tp - tp0;
+ end;
+ end
+ n_nnz = length(sum(is_constraint,2));
+ if n_nnz > 0
+ constraint_index = cell(n_nnz,1);
+ for i= 1:n_nnz
+ constraint_index{i} = find(is_constraint(i,:));
+ end;
+ end;
+ if shocks_present
+ n_periods = length(shock_periods);
+ shock_paths_cell = shock_paths;
+ clear shock_paths;
+ shock_paths = zeros(n_periods, size(total_periods));
+ is_shock = zeros(size(total_periods), n_periods);
+ for i = 1:n_periods
+ period_i = shock_periods{i};
+ %period_i
+ tp = total_periods(1);
+ if size(period_i) > 1
+ init_periods = period_i(1);
+ tp_end = period_i(end);
+ else
+ init_periods = period_i;
+ tp_end = period_i;
+ end;
+ tp0 = tp;
+ while tp < init_periods
+ tp = tp + 1;
+ end
+ j = 0;
+ while tp <= tp_end
+ is_shock(tp - tp0 + 1, i) = 1;
+ shock_paths(i, tp - tp0 + 1) = shock_paths_cell{i}(j + 1);
+ tp = tp + 1;
+ j = j + 1;
+ end;
+ if tp - tp0 > max_periods_simulation
+ max_periods_simulation = tp - tp0;
+ end;
+ end;
+ n_nnz = length(sum(is_shock,2));
+ if n_nnz > 0
+ shock_index = cell(n_nnz, 1);
+ for i= 1:n_nnz
+ shock_index{i} = find(is_shock(i,:));
+ end;
+ end
+ end
+else
+ is_constraint = ones(size(constrained_paths));
+end
+
+
+
+maximum_lag = M_.maximum_lag;
+
+ys = oo_.steady_state;
+ny = size(ys,1);
+xs = [oo_.exo_steady_state ; oo_.exo_det_steady_state];
+nx = size(xs,1);
+
+constrained_periods = max_periods_simulation;
+n_endo_constrained = size(constrained_vars,1);
+if isfield(options_cond_fcst,'controlled_varexo')
+ n_control_exo = size(options_cond_fcst.controlled_varexo, 1);
+ if n_control_exo ~= n_endo_constrained
+ error(['det_cond_forecast: the number of exogenous controlled variables (' int2str(n_control_exo) ') has to be equal to the number of constrained endogenous variabes (' int2str(n_endo_constrained) ')'])
+ end;
+else
+ error('det_cond_forecast: to run a deterministic conditional forecast you have to specified the exogenous variables controlled using the option controlled_varexo in forecast command');
+end;
+
+% if n_endo_constrained == 0
+% options_.ep.use_bytecode = options_.bytecode;
+% data_set = extended_path(initial_conditions, max_periods_simulation);
+% end
+
+if length(varargin) >= 1
+ controlled_varexo = options_cond_fcst.controlled_varexo;
+else
+ exo_names = M_.exo_names;
+ controlled_varexo = zeros(1,n_control_exo);
+ for i = 1:nx
+ for j=1:n_control_exo
+ if strcmp(deblank(exo_names(i,:)), deblank(options_cond_fcst.controlled_varexo(j,:)))
+ controlled_varexo(j) = i;
+ end
+ end
+ end
+end
+
+%todo check if zero => error message
+
+save_options_initval_file = options_.initval_file;
+options_.initval_file = '__';
+
+[pos_constrained_pf, junk] = find(constrained_perfect_foresight);
+indx_endo_solve_pf = constrained_vars(pos_constrained_pf);
+if isempty(indx_endo_solve_pf)
+ pf = 0;
+else
+ pf = length(indx_endo_solve_pf);
+end;
+indx_endo_solve_surprise = setdiff(constrained_vars, indx_endo_solve_pf);
+
+if isempty(indx_endo_solve_surprise)
+ surprise = 0;
+else
+ surprise = length(indx_endo_solve_surprise);
+end;
+
+
+eps = options_.solve_tolf;
+maxit = options_.simul.maxit;
+
+
+past_val = 0;
+save_options_periods = options_.periods;
+options_.periods = options_cond_fcst.periods;
+save_options_dynatol_f = options_.dynatol.f;
+options_.dynatol.f = 1e-8;
+eps1 = 1e-7;%1e-4;
+exo = zeros(maximum_lag + options_cond_fcst.periods, nx);
+endo = zeros(maximum_lag + options_cond_fcst.periods, ny);
+endo(1,:) = oo_.steady_state';
+verbosity = options_.verbosity;
+options_.verbosity = 0;
+
+
+% if all the endogenous paths are perfectly anticipated we do not need to
+% implement the extended path
+if pf && ~surprise
+ time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods;
+ second_system_size = pf * constrained_periods;
+ J = zeros(second_system_size,second_system_size);
+ r = zeros(second_system_size,1);
+ indx_endo = zeros(second_system_size,1);
+ col_count = 1;
+ for j = 1:length(constrained_vars)
+ indx_endo(col_count : col_count + constrained_periods - 1) = constrained_vars(j) + (time_index_constraint - 1) * ny;
+ col_count = col_count + constrained_periods;
+ end;
+ make_ex_;
+ make_y_;
+ it = 1;
+ convg = 0;
+ normra = 1e+50;
+ while ~convg && it <= maxit
+ disp('---------------------------------------------------------------------------------------------');
+ disp(['iteration ' int2str(it)]);
+ not_achieved = 1;
+ alpha = 1;
+ while not_achieved
+ simul();
+ result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * constrained_periods;
+ if result
+ y = oo_.endo_simul(constrained_vars, time_index_constraint);
+ ys = oo_.endo_simul(indx_endo);
+ col_count = 1;
+ for j = 1:length(constrained_vars)
+ y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
+ r(col_count:col_count+constrained_periods - 1) = (y - constrained_paths(j,1:constrained_periods))';
+ col_count = col_count + constrained_periods;
+ end;
+ normr = norm(r, 1);
+ end;
+ if (~oo_.deterministic_simulation.status || ~result || normr > normra) && it > 1
+ not_achieved = 1;
+ alpha = alpha / 2;
+ col_count = 1;
+ for j = controlled_varexo'
+ oo_.exo_simul(time_index_constraint,j) = (old_exo(:,j) + alpha * D_exo(col_count: (col_count + constrained_periods - 1)));
+ col_count = col_count + constrained_periods;
+ end;
+ disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha)]);
+ oo_.endo_simul = repmat(oo_.steady_state, 1, options_cond_fcst.periods + 2);
+ else
+ not_achieved = 0;
+ end;
+ end;
+
+ per = 30;
+ z = oo_.endo_simul(:, 1 : per + 2 );
+ zx = oo_.exo_simul(1: per + 2,:);
+ g1 = spalloc(M_.endo_nbr * (per ), M_.endo_nbr * (per ), 3* M_.endo_nbr * per );
+ g1_x = spalloc(M_.endo_nbr * (per ), M_.exo_nbr, M_.endo_nbr * (per )* M_.exo_nbr );
+ lag_indx = find(M_.lead_lag_incidence(1,:));
+ cur_indx = M_.endo_nbr + find(M_.lead_lag_incidence(2,:));
+ lead_indx = 2 * M_.endo_nbr + find(M_.lead_lag_incidence(3,:));
+ cum_l1 = 0;
+ cum_index_d_y_x = [];
+ indx_x = [];
+ for k = 1 : per
+ if k == 1
+ if (isfield(M_,'block_structure'))
+ data1 = M_.block_structure.block;
+ Size = length(M_.block_structure.block);
+ else
+ data1 = M_;
+ Size = 1;
+ end;
+ data1 = M_;
+ if (options_.bytecode)
+ [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
+ else
+ [zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, oo_.steady_state, k);
+ data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
+ data1.g1 = g1b(:,1 : end - M_.exo_nbr);
+ chck = 0;
+ end;
+ mexErrCheck('bytecode', chck);
+ end;
+ if k == 1
+ g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end);
+ elseif k == per
+ g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k,M_.endo_nbr * (k -2) + [lag_indx cur_indx]) = data1.g1(:,1:M_.nspred + M_.endo_nbr);
+ else
+ g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k, M_.endo_nbr * (k -2) + [lag_indx cur_indx lead_indx]) = data1.g1;
+ end;
+ l2 = 1;
+ pf_c = 1;
+ if k <= constrained_periods
+ for l = constraint_index{k}
+ l1 = controlled_varexo(l);
+ g1_x(M_.endo_nbr * (k - 1) + 1:M_.endo_nbr * k,1 + cum_l1) = data1.g1_x(:,l1);
+ if k == 1
+ indx_x(l2) = l ;
+ l2 = l2 + 1;
+ for ii = 2:constrained_periods
+ indx_x(l2) = length(controlled_varexo) + pf * (ii - 2) + constraint_index{k + ii - 1}(pf_c);
+ l2 = l2 + 1;
+ end;
+ pf_c = pf_c + 1;
+ cum_index_d_y_x = [cum_index_d_y_x; constrained_vars(l)];
+ else
+ cum_index_d_y_x = [cum_index_d_y_x; constrained_vars(l) + (k - 1) * M_.endo_nbr];
+ end
+ cum_l1 = cum_l1 + length(l1);
+ end;
+ end;
+ end;
+
+ d_y_x = - g1 \ g1_x;
+
+ cum_l1 = 0;
+ count_col = 1;
+ cum_index_J = 1:length(cum_index_d_y_x(indx_x));
+ J= zeros(length(cum_index_J));
+ for j1 = 1:length(controlled_varexo)
+ cum_l1 = 0;
+ for k = 1:(constrained_periods)
+ l1 = constraint_index{k};
+ l1 = find(constrained_perfect_foresight(l1) | (k == 1));
+ if constraint_index{k}( j1)
+ J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col));
+ count_col = count_col + 1;
+ end
+ cum_l1 = cum_l1 + length(l1);
+ end
+ cum_l1 = cum_l1 + length(constrained_vars(j1));
+ end;
+
+
+% col_count = 1;
+% for j = controlled_varexo'
+% for time = time_index_constraint
+% saved = oo_.exo_simul(time,j);
+% oo_.exo_simul(time,j) = oo_.exo_simul(time,j) + eps1;
+% simul();
+% J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
+% oo_.exo_simul(time,j) = saved;
+% col_count = col_count + 1;
+% end;
+% end;
+% J1
+% sdfmlksdf;
+
+ disp(['iteration ' int2str(it) ' error = ' num2str(normr)]);
+
+ if normr <= eps
+ convg = 1;
+ disp('convergence achieved');
+ else
+ % Newton update on exogenous shocks
+ old_exo = oo_.exo_simul(time_index_constraint,:);
+ D_exo = - J \ r;
+ col_count = 1;
+ %constrained_periods
+ for j = controlled_varexo'
+ oo_.exo_simul(time_index_constraint,j) = oo_.exo_simul(time_index_constraint,j) + D_exo(col_count: (col_count + constrained_periods - 1));
+ col_count = col_count + constrained_periods - 1;
+ end;
+ end;
+ it = it + 1;
+ normra = normr;
+ end;
+ endo = oo_.endo_simul';
+ exo = oo_.exo_simul;
+else
+ for t = 1:constrained_periods
+
+ if direct_mode && ~isempty(is_constraint)
+ [pos_constrained_pf, junk] = find(constrained_perfect_foresight .* is_constraint(t, :)');
+ indx_endo_solve_pf = constrained_vars(pos_constrained_pf);
+ if isempty(indx_endo_solve_pf)
+ pf = 0;
+ else
+ pf = length(indx_endo_solve_pf);
+ end;
+
+ [pos_constrained_surprise, junk] = find((1-constrained_perfect_foresight) .* is_constraint(t, :)');
+ indx_endo_solve_surprise = constrained_vars(pos_constrained_surprise);
+
+ if isempty(indx_endo_solve_surprise)
+ surprise = 0;
+ else
+ surprise = length(indx_endo_solve_surprise);
+ end;
+ end;
+
+ if direct_mode && ~isempty(is_shock)
+ [pos_shock_pf, junk] = find(shock_perfect_foresight .* is_shock(t, :)');
+ indx_endo_solve_pf = shock_vars(pos_shock_pf);
+ if isempty(indx_endo_solve_pf)
+ b_pf = 0;
+ else
+ b_pf = length(indx_endo_solve_pf);
+ end;
+
+ [pos_shock_surprise, junk] = find((1-shock_perfect_foresight) .* is_shock(t, :)');
+ indx_endo_solve_surprise = shock_vars(pos_shock_surprise);
+
+ if isempty(indx_endo_solve_surprise)
+ b_surprise = 0;
+ else
+ b_surprise = length(indx_endo_solve_surprise);
+ end;
+ end;
+
+ disp('===============================================================================================');
+ disp(['t=' int2str(t) ' conditional (surprise=' int2str(surprise) ' perfect foresight=' int2str(pf) ') unconditional (surprise=' int2str(b_surprise) ' perfect foresight=' int2str(b_pf) ')']);
+ disp('===============================================================================================');
+ if t == 1
+ make_ex_;
+ if maximum_lag > 0
+ exo_init = oo_.exo_simul;
+ else
+ exo_init = zeros(size(oo_.exo_simul));
+ end
+ make_y_;
+ end;
+ %exo_init
+ oo_.exo_simul = exo_init;
+ oo_.endo_simul(:,1) = initial_conditions;
+
+ time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods - t + 1;
+
+ if direct_mode
+ nb_shocks = length(plan.shock_vars_);
+ if nb_shocks > 0
+ shock_index_t = shock_index{t};
+ shock_vars_t = shock_vars(shock_index_t);
+ for shock_indx = shock_index_t
+ if shock_perfect_foresight(shock_indx)
+ oo_.exo_simul(time_index_constraint,shock_vars_t(shock_indx)) = shock_paths(shock_indx,t:constrained_periods);
+ else
+ oo_.exo_simul(maximum_lag + 1,shock_vars_t(shock_indx)) = shock_paths(shock_indx,t);
+ end
+ end
+ end
+ end
+ %Compute the initial path using the first order solution around the
+ % steady-state
+ oo_.endo_simul(:, 1) = oo_.endo_simul(:, 1) - oo_.steady_state;
+ for jj = 2 : (options_.periods + 2)
+ oo_.endo_simul(:,jj) = oo_.dr.ghx(oo_.dr.inv_order_var,:) * oo_.endo_simul(oo_.dr.state_var, jj-1) + oo_.dr.ghu * oo_.exo_simul(jj,:)';
+ end
+ for jj = 1 : (options_.periods + 2)
+ oo_.endo_simul(:,jj) = oo_.steady_state + oo_.endo_simul(:,jj);
+ end
+
+ constraint_index_t = constraint_index{t};
+ controlled_varexo_t = controlled_varexo(constraint_index_t);
+ constrained_vars_t = constrained_vars(constraint_index_t);
+
+ second_system_size = surprise + pf * (constrained_periods - t + 1);
+ indx_endo = zeros(second_system_size,1);
+ col_count = 1;
+ for j = 1:length(constrained_vars_t)
+ if constrained_perfect_foresight(j)
+ indx_endo(col_count : col_count + constrained_periods - t) = constrained_vars(j) + (time_index_constraint - 1) * ny;
+ col_count = col_count + constrained_periods - t + 1;
+ else
+ indx_endo(col_count) = constrained_vars(j) + maximum_lag * ny;
+ col_count = col_count + 1;
+ end;
+ end;
+
+ r = zeros(second_system_size,1);
+
+ convg = 0;
+ it = 1;
+ while ~convg && it <= maxit
+ disp('-----------------------------------------------------------------------------------------------');
+ disp(['iteration ' int2str(it) ' time ' int2str(t)]);
+ not_achieved = 1;
+ alpha = 1;
+ while not_achieved
+ simul();
+ result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * length(time_index_constraint);
+ if (~oo_.deterministic_simulation.status || ~result) && it > 1
+ not_achieved = 1;
+ alpha = alpha / 2;
+ col_count = 1;
+ for j1 = constraint_index_t
+ j = controlled_varexo(j1);
+ if constrained_perfect_foresight(j1)
+ oo_.exo_simul(time_index_constraint,j) = (old_exo(time_index_constraint,j) + alpha * D_exo(col_count: col_count + constrained_periods - t));
+ col_count = col_count + constrained_periods - t + 1;
+ else
+ oo_.exo_simul(maximum_lag + 1,j) = old_exo(maximum_lag + 1,j) + alpha * D_exo(col_count);
+ col_count = col_count + 1;
+ end;
+ end;
+ disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha) ' result=' num2str(result) ' sum(sum)=' num2str(sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint))))) ' ny * length(time_index_constraint)=' num2str(ny * length(time_index_constraint)) ' oo_.deterministic_simulation.status=' num2str(oo_.deterministic_simulation.status)]);
+ oo_.endo_simul = [initial_conditions repmat(oo_.steady_state, 1, options_cond_fcst.periods + 1)];
+ else
+ not_achieved = 0;
+ end;
+ end;
+ if t==constrained_periods
+ ycc = oo_.endo_simul;
+ end;
+ yc = oo_.endo_simul(:,maximum_lag + 1);
+ ys = oo_.endo_simul(indx_endo);
+
+ col_count = 1;
+ for j = constraint_index_t
+ if constrained_perfect_foresight(j)
+ y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
+ r(col_count:col_count+constrained_periods - t) = (y - constrained_paths(j,t:constrained_periods))';
+ col_count = col_count + constrained_periods - t + 1;
+ else
+ y = yc(constrained_vars(j));
+ r(col_count) = y - constrained_paths(j,t);
+ col_count = col_count + 1;
+ end;
+ end
+
+ disp('computation of derivatives w.r. to exogenous shocks');
+
+ per = 30;
+ z = oo_.endo_simul(:, 1 : per + 2 );
+ zx = oo_.exo_simul(1: per + 2,:);
+ g1 = spalloc(M_.endo_nbr * (per ), M_.endo_nbr * (per ), 3* M_.endo_nbr * per );
+ g1_x = spalloc(M_.endo_nbr * (per ), M_.exo_nbr, M_.endo_nbr * (per )* M_.exo_nbr );
+ lag_indx = find(M_.lead_lag_incidence(1,:));
+ cur_indx = M_.endo_nbr + find(M_.lead_lag_incidence(2,:));
+ lead_indx = 2 * M_.endo_nbr + find(M_.lead_lag_incidence(3,:));
+ cum_l1 = 0;
+ %indx_x = zeros(length(constraint_index_t)* constrained_periods, 1);
+ cum_index_d_y_x = [];
+ indx_x = [];
+ for k = 1 : per
+ if k == 1
+ if (isfield(M_,'block_structure'))
+ data1 = M_.block_structure.block;
+ Size = length(M_.block_structure.block);
+ else
+ data1 = M_;
+ Size = 1;
+ end;
+ data1 = M_;
+ if (options_.bytecode)
+ [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
+ else
+ [zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, oo_.steady_state, k);
+ data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
+ data1.g1 = g1b(:,1 : end - M_.exo_nbr);
+ chck = 0;
+ end;
+ mexErrCheck('bytecode', chck);
+ end;
+ if k == 1
+ g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end);
+ elseif k == per
+ g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k,M_.endo_nbr * (k -2) + [lag_indx cur_indx]) = data1.g1(:,1:M_.nspred + M_.endo_nbr);
+ else
+ g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k, M_.endo_nbr * (k -2) + [lag_indx cur_indx lead_indx]) = data1.g1;
+ end;
+
+ l2 = 1;
+ pf_c = 1;
+ if t + k - 1 <= constrained_periods
+ for l = constraint_index{t + k - 1}
+ l1 = controlled_varexo(l);
+ if (k == 1) || ((k > 1) && constrained_perfect_foresight(l))
+ g1_x(M_.endo_nbr * (k - 1) + 1:M_.endo_nbr * k,1 + cum_l1) = data1.g1_x(:,l1);
+ if k == 1
+ if constrained_perfect_foresight(l)
+ indx_x(l2) = l ;
+ l2 = l2 + 1;
+ for ii = 2:constrained_periods - t + 1
+ indx_x(l2) = length(controlled_varexo) + pf * (ii - 2) + constraint_index{k + ii - 1}(pf_c);
+ l2 = l2 + 1;
+ end;
+ pf_c = pf_c + 1;
+ else
+ indx_x(l2) = l;
+ l2 = l2 + 1;
+ end;
+ cum_index_d_y_x = [cum_index_d_y_x; constrained_vars_t(l)];
+ else
+ cum_index_d_y_x = [cum_index_d_y_x; constrained_vars_t(l) + (k - 1) * M_.endo_nbr];
+ end
+ cum_l1 = cum_l1 + length(l1);
+ end;
+ end;
+ end;
+ end;
+
+ d_y_x = - g1 \ g1_x;
+
+
+ cum_l1 = 0;
+ count_col = 1;
+ cum_index_J = 1:length(cum_index_d_y_x(indx_x));
+ J= zeros(length(cum_index_J));
+ for j1 = constraint_index_t
+ if constrained_perfect_foresight(j1)
+ cum_l1 = 0;
+ for k = 1:(constrained_periods - t + 1)
+ l1 = constraint_index{k};
+ l1 = find(constrained_perfect_foresight(l1) | (k == 1));
+ if constraint_index{k}( j1)
+ J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col));
+ count_col = count_col + 1;
+ end
+ cum_l1 = cum_l1 + length(l1);
+ end
+ else
+ J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col));
+ count_col = count_col + 1;
+ end
+ cum_l1 = cum_l1 + length(constrained_vars_t(j1));
+ end;
+
+
+% % Numerical computation of the derivatives in the second systme
+% J1 = [];
+% col_count = 1;
+% for j = constraint_index_t
+% j_pos = controlled_varexo(j);
+% if constrained_perfect_foresight(j)
+% for time = time_index_constraint
+% saved = oo_.exo_simul(time,j_pos);
+% oo_.exo_simul(time,j_pos) = oo_.exo_simul(time,j_pos) + eps1;
+% simul();
+% J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
+% oo_.exo_simul(time,j_pos) = saved;
+% col_count = col_count + 1;
+% end;
+% else
+% saved = oo_.exo_simul(maximum_lag+1,j_pos);
+% oo_.exo_simul(maximum_lag+1,j_pos) = oo_.exo_simul(maximum_lag+1,j_pos) + eps1;
+% simul();
+% % indx_endo
+% J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
+% % J(:,col_count) = (oo_.endo_simul((pp - 1) * M_.endo_nbr + 1: pp * M_.endo_nbr) - ys) / eps1;
+% oo_.exo_simul(maximum_lag+1,j_pos) = saved;
+% col_count = col_count + 1;
+% end;
+% end;
+% disp('J1');
+% disp(full(J1));
+% sdfmlk;
+
+
+ normr = norm(r, 1);
+
+ disp(['iteration ' int2str(it) ' error = ' num2str(normr) ' at time ' int2str(t)]);
+
+ if normr <= eps
+ convg = 1;
+ disp('convergence achieved');
+ else
+ % Newton update on exogenous shocks
+ try
+ D_exo = - J \ r;
+ catch
+ [V, D] = eig(full(J));
+ ev = diag(D);
+ [ev abs(ev)]
+ z_root = find(abs(ev) < 1e-4);
+ z_root
+ disp(V(:,z_root));
+ end;
+ old_exo = oo_.exo_simul;
+ col_count = 1;
+ for j = constraint_index_t
+ j_pos=controlled_varexo(j);
+ if constrained_perfect_foresight(j)
+ oo_.exo_simul(time_index_constraint,j_pos) = (oo_.exo_simul(time_index_constraint,j_pos) + D_exo(col_count:(col_count + constrained_periods - t) ));
+ col_count = col_count + constrained_periods - t + 1;
+ else
+ oo_.exo_simul(maximum_lag + 1,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos) + D_exo(col_count);
+ col_count = col_count + 1;
+ end;
+ end;
+ end;
+ it = it + 1;
+ end;
+ if ~convg
+ error(['convergence not achived at time ' int2str(t) ' after ' int2str(it) ' iterations']);
+ end;
+ for j = constraint_index_t
+ j_pos=controlled_varexo(j);
+ if constrained_perfect_foresight(j)
+ % in case of mixed surprise and perfect foresight on the
+ % endogenous path, at each date all the exogenous paths have to be
+ % stored. The paths are stacked in exo.
+ for time = time_index_constraint;
+ exo(past_val + time,j_pos) = oo_.exo_simul(time,j_pos);
+ end
+ else
+ exo(maximum_lag + t,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos);
+ end;
+ end;
+ past_val = past_val + length(time_index_constraint);
+ if t < constrained_periods
+ endo(maximum_lag + t,:) = yc;
+ else
+ endo(maximum_lag + t :maximum_lag + options_cond_fcst.periods ,:) = ycc(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods - constrained_periods + 1)';
+ end;
+ initial_conditions = yc;
+ if maximum_lag > 0
+ exo_init(1,:) = exo(maximum_lag + t,:);
+ end;
+ end;
+end;
+options_.periods = save_options_periods;
+options_.dynatol.f = save_options_dynatol_f;
+options_.initval_file = save_options_initval_file;
+options_.verbosity = verbosity;
+oo_.endo_simul = endo';
+oo_.exo_simul = exo;
+if direct_mode
+ data_set = dseries([endo(2:end,1:M_.orig_endo_nbr) exo(2:end,:)], total_periods(1), {plan.endo_names{:} plan.exo_names{:}}, {plan.endo_names{:} plan.exo_names{:}});
+end;
diff --git a/matlab/flip_plan.m b/matlab/flip_plan.m
new file mode 100644
index 0000000000000000000000000000000000000000..6859e76a8504db362529cbf4fea8aaa1a79ef22a
--- /dev/null
+++ b/matlab/flip_plan.m
@@ -0,0 +1,89 @@
+function plan = flip_plan(plan, exogenous, endogenous, expectation_type, date, value)
+% Adds to the forecast scenario a conditional forecast shock (the path of an endogenous variable is constrained and the values compatible values of the related exogenous variable will be compued)
+%
+% INPUTS
+% o plan [structure] A structure describing the different shocks and the implied variables, the date of the shocks and the path of the shock (forecast scenario).
+% The plan structure is created by the functions init_plan, basic_plan and flip_plan
+% o exogenous [string] A string containing the name of the endogenous variable with a constrained path.
+% o endogenous [string] A string containing the name of the exogenous. This exogenous variable will be en endogenous variable when the conditional forecast will be perform.
+% o expectation_type [string] A string indicating the type of expectation: 'surprise' for an unexpected shock, and 'perfect_foresight' for a perfectly anticpated shock.
+% o date [dates] The period of the shock
+% o value [array of double] A vector of double containing the constrined path on the endogenous variable
+%
+%
+% OUTPUTS
+% plan [structure] Returns a structure containing the updated forecast scenario.
+%
+%
+% Copyright (C) 2013-2014 Dynare Team
+%
+% 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/>.
+ if ~ischar(expectation_type) || size(expectation_type,1) ~= 1
+ error(['in flip_plan the fourth argument should be a string containing the simulation type (''perfect_foresight'' or ''surprise'')']);
+ end
+ exogenous = strtrim(exogenous);
+ ix = find(strcmp(exogenous, plan.endo_names));
+ if isempty(ix)
+ error(['in flip_plan the second argument ' exogenous ' is not an endogenous variable']);
+ end;
+ endogenous = strtrim(endogenous);
+ iy = find(strcmp(endogenous, plan.exo_names));
+ if isempty(iy)
+ error(['in flip_plan the third argument ' endogenous ' is not an exogenous variable']);
+ end;
+ sdate = length(date);
+ if sdate > 1
+ if date(1) < plan.date(1) || date(end) > plan.date(end)
+ error(['in flip_plan the fifth argument (date=' date ') must lay inside the plan.date ' plan.date]);
+ end
+ else
+ if date < plan.date(1) || date > plan.date(end)
+ error(['in flip_plan the fifth argument (date=' date ') must lay iside the plan.date ' plan.date]);
+ end
+ end
+ if ~isempty(plan.shock_vars_)
+ common_var = find(iy == plan.shock_vars_);
+ if ~isempty(common_var)
+ common_date = intersect(date, plan.shock_date_{common_var});
+ if ~isempty(common_date)
+ if common_date.length > 1
+ the_dates = [cell2mat(strings(common_date(1))) ':' cell2mat(strings(common_date(end)))];
+ else
+ the_dates = cell2mat(strings(common_date));
+ end
+ error(['Impossible case: ' plan.exo_names{plan.shock_vars_(common_var)} ' is used both as a shock and as an endogenous variable to control the path of ' plan.endo_names{ix} ' at the dates ' the_dates]);
+ end
+ end
+ end
+ if isempty(plan.constrained_vars_)
+ plan.constrained_vars_ = ix;
+ plan.options_cond_fcst_.controlled_varexo = iy;
+ if strcmp(expectation_type, 'perfect_foresight')
+ plan.constrained_perfect_foresight_ = 1;
+ else
+ plan.constrained_perfect_foresight_ = 0;
+ end
+ else
+ plan.constrained_vars_ = [plan.constrained_vars_ ; ix];
+ plan.options_cond_fcst_.controlled_varexo = [plan.options_cond_fcst_.controlled_varexo ; iy];
+ if strcmp(expectation_type, 'perfect_foresight')
+ plan.constrained_perfect_foresight_ = [plan.constrained_perfect_foresight_ ; 1];
+ else
+ plan.constrained_perfect_foresight_ = [plan.constrained_perfect_foresight_ ; 0];
+ end
+ end
+ plan.constrained_date_{length(plan.constrained_date_) + 1} = date;
+ plan.constrained_paths_{length(plan.constrained_paths_) + 1} = value;
diff --git a/matlab/init_plan.m b/matlab/init_plan.m
new file mode 100644
index 0000000000000000000000000000000000000000..45be92d460b14a34c9bce5da91f724493320fc2d
--- /dev/null
+++ b/matlab/init_plan.m
@@ -0,0 +1,47 @@
+function plan = init_plan(date)
+% Creates and initializes a new forecast scenario
+%
+% INPUTS
+% o date [dates] The period of the forecast
+%
+%
+% OUTPUTS
+% plan [structure] Returns a structure containing a new forecast scenario
+%
+%
+% Copyright (C) 2013-2014 Dynare Team
+%
+% 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/>.
+ global M_
+ plan = struct();
+ plan.date = date;
+ endo_names_length = size(M_.endo_names,2);
+ plan.endo_names = deblank(mat2cell(M_.endo_names(1:M_.orig_endo_nbr,:),ones(1,M_.orig_endo_nbr),endo_names_length));
+ exo_names_length = size(M_.exo_names,2);
+ plan.exo_names = deblank(mat2cell(M_.exo_names(1:M_.exo_nbr,:),ones(1,M_.exo_nbr),exo_names_length));
+ plan.constrained_vars_ = [];
+ plan.constrained_paths_ = [];
+ plan.constrained_date_ = [];
+ plan.constrained_perfect_foresight_ = [];
+ plan.shock_vars_ = [];
+ plan.shock_paths_ = [];
+ plan.shock_date_ = [];
+ plan.shock_perfect_foresight_ = [];
+ plan.options_cond_fcst_ = struct();
+ plan.options_cond_fcst_.parameter_set = 'calibration';
+ plan.options_cond_fcst_.simulation_type = 'deterministic';
+ plan.options_cond_fcst_.controlled_varexo = [];
+
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index b6ec34311894ec8139a39451a5d92b21c0a9ce1f..c04819b7213f5f9e6a8ceda4d2485e1ad717b14b 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -89,7 +89,7 @@ class ParsingDriver;
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
%token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED
%token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
-%token DATAFILE FILE DETERMINISTIC DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
+%token DATAFILE FILE DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
%token <string_val> FLOAT_NUMBER DATES
@@ -115,8 +115,8 @@ class ParsingDriver;
%token <string_val> QUOTED_STRING
%token QZ_CRITERIUM QZ_ZERO_THRESHOLD FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE
%token RELATIVE_IRF REPLIC SIMUL_REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY
-%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SIMULATION_TYPE ENDOGENOUS_TERMINAL_PERIOD
-%token SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL STOCHASTIC SOLVE_ALGO SOLVER_PERIODS
+%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED ENDOGENOUS_TERMINAL_PERIOD
+%token SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO SOLVER_PERIODS
%token STDERR STEADY STOCH_SIMUL SURPRISE SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME
%token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY DISCRETIONARY_TOL
%token <string_val> TEX_NAME
@@ -877,12 +877,6 @@ period_list : period_list COMMA INT_NUMBER
{ driver.add_period($1); }
;
-expectation_type : PERFECT_FORESIGHT
- { driver.add_expectation_pf(true); }
- | SURPRISE
- { driver.add_expectation_pf(false); }
- ;
-
sigma_e : SIGMA_E EQUAL '[' triangular_matrix ']' ';' { driver.do_sigma_e(); };
value_list : value_list COMMA '(' expression ')'
@@ -2233,7 +2227,6 @@ conditional_forecast_option : o_periods
| o_conf_sig
| o_controlled_varexo
| o_parameter_set
- | o_simulation_type
;
plot_conditional_forecast : PLOT_CONDITIONAL_FORECAST symbol_list ';'
@@ -2252,8 +2245,6 @@ conditional_forecast_paths_shock_list : conditional_forecast_paths_shock_elem
conditional_forecast_paths_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
{ driver.add_det_shock($2, true); }
- | VAR symbol ';' PERIODS period_list ';' VALUES value_list ';' EXPECTATION expectation_type ';'
- { driver.add_det_shock($2, true); }
;
steady_state_model : STEADY_STATE_MODEL ';' { driver.begin_steady_state_model(); }
@@ -2563,11 +2554,6 @@ o_parameter_set : PARAMETER_SET EQUAL PRIOR_MODE
| PARAMETER_SET EQUAL CALIBRATION
{ driver.option_str("parameter_set", "calibration"); }
;
-o_simulation_type : SIMULATION_TYPE EQUAL DETERMINISTIC
- { driver.option_str("simulation_type", "deterministic"); }
- | SIMULATION_TYPE EQUAL STOCHASTIC
- { driver.option_str("simulation_type", "stochastic"); }
- ;
o_ms_drop : DROP EQUAL INT_NUMBER { driver.option_num("ms.drop", $3); };
o_ms_mh_replic : MH_REPLIC EQUAL INT_NUMBER { driver.option_num("ms.mh_replic", $3); };
o_freq : FREQ EQUAL INT_NUMBER
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 98dda2663c4dda974eeda8e15baebe0113f45b41..f8395e20cdb6c547520c7fb973229fc9fd52fb78 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -505,9 +505,6 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>posterior_mode {return token::POSTERIOR_MODE; }
<DYNARE_STATEMENT>posterior_mean {return token::POSTERIOR_MEAN; }
<DYNARE_STATEMENT>posterior_median {return token::POSTERIOR_MEDIAN; }
-<DYNARE_STATEMENT>simulation_type {return token::SIMULATION_TYPE; }
-<DYNARE_STATEMENT>deterministic {return token::DETERMINISTIC; }
-<DYNARE_STATEMENT>stochastic {return token::STOCHASTIC; }
<DYNARE_STATEMENT>k_order_solver {return token::K_ORDER_SOLVER; }
<DYNARE_STATEMENT>filter_covariance {return token::FILTER_COVARIANCE; }
<DYNARE_STATEMENT>filter_decomposition {return token::FILTER_DECOMPOSITION; }
@@ -537,8 +534,6 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_BLOCK>stderr {return token::STDERR;}
<DYNARE_BLOCK>values {return token::VALUES;}
<DYNARE_BLOCK>corr {return token::CORR;}
-<DYNARE_BLOCK>surprise {return token::SURPRISE;}
-<DYNARE_BLOCK>perfect_foresight {return token::PERFECT_FORESIGHT;}
<DYNARE_BLOCK>periods {return token::PERIODS;}
<DYNARE_BLOCK>cutoff {return token::CUTOFF;}
<DYNARE_BLOCK>mfs {return token::MFS;}
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index c556220f03684f29cda21a15e8d42b8601a66787..f4bf3a755226e3cbbc158e59f96146e20231b7cf 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -670,12 +670,10 @@ ParsingDriver::add_det_shock(string *var, bool conditional_forecast)
v.push_back(dse);
}
- det_shocks[symb_id].first = v;
- det_shocks[symb_id].second = det_shocks_expectation_pf;
-
+ det_shocks[symb_id] = v;
+
det_shocks_periods.clear();
det_shocks_values.clear();
- det_shocks_expectation_pf = false;
delete var;
}
@@ -795,12 +793,6 @@ ParsingDriver::add_value(string *v)
det_shocks_values.push_back(id);
}
-void
-ParsingDriver::add_expectation_pf(bool pf)
-{
- det_shocks_expectation_pf = pf;
-}
-
void
ParsingDriver::begin_svar_identification()
{
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index 017c9361061d5264ae8080064da2ada104c351fa..e856a2e0fe24e80748fca68f3d8a4215196c27b9 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -134,8 +134,6 @@ private:
vector<pair<int, int> > det_shocks_periods;
//! Temporary storage for values of deterministic shocks
vector<expr_t> det_shocks_values;
- //! Temporary storage for perfect foresight of deterministic shocks in conditional forecast
- bool det_shocks_expectation_pf;
//! Temporary storage for variances of shocks
ShocksStatement::var_and_std_shocks_t var_shocks;
//! Temporary storage for standard errors of shocks
@@ -215,7 +213,7 @@ private:
bool nostrict;
public:
- ParsingDriver(WarningConsolidation &warnings_arg, bool nostrict_arg) : det_shocks_expectation_pf(false), warnings(warnings_arg), nostrict(nostrict_arg) { };
+ ParsingDriver(WarningConsolidation &warnings_arg, bool nostrict_arg) : warnings(warnings_arg), nostrict(nostrict_arg) { };
//! Starts parsing, and constructs the MOD file representation
/*! The returned pointer should be deleted after use */
@@ -347,8 +345,6 @@ public:
//! Adds a deterministic shock value
/*! \param v a string containing a (possibly negative) numeric constant */
void add_value(string *v);
- //! Adds a expectation type for conditional forecast with deterministic simulation
- void add_expectation_pf(bool pf);
//! Writes a Sigma_e block
void do_sigma_e();
//! Ends row of Sigma_e block
diff --git a/preprocessor/Shocks.cc b/preprocessor/Shocks.cc
index 4b4037b2f8f32cfd62dac66aea4d15dbb7db317e..73032fbb540e059ecf48a108960c3023ad16a9dd 100644
--- a/preprocessor/Shocks.cc
+++ b/preprocessor/Shocks.cc
@@ -44,11 +44,11 @@ AbstractShocksStatement::writeDetShocks(ostream &output) const
bool exo_det = (symbol_table.getType(it->first) == eExogenousDet);
int set_shocks_index = ((int) mshocks) + 2 * ((int) exo_det);
- for (size_t i = 0; i < it->second.first.size(); i++)
+ for (size_t i = 0; i < it->second.size(); i++)
{
- const int &period1 = it->second.first[i].period1;
- const int &period2 = it->second.first[i].period2;
- const expr_t value = it->second.first[i].value;
+ const int &period1 = it->second[i].period1;
+ const int &period2 = it->second[i].period2;
+ const expr_t value = it->second[i].value;
if (period1 == period2)
{
@@ -339,7 +339,7 @@ ConditionalForecastPathsStatement::checkPass(ModFileStructure &mod_file_struct,
it != paths.end(); it++)
{
int this_path_length = 0;
- const vector<AbstractShocksStatement::DetShockElement> &elems = it->second.first;
+ const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
for (int i = 0; i < (int) elems.size(); i++)
// Period1 < Period2, as enforced in ParsingDriver::add_period()
this_path_length = max(this_path_length, elems[i].period2);
@@ -368,16 +368,14 @@ ConditionalForecastPathsStatement::writeOutput(ostream &output, const string &ba
if (it == paths.begin())
{
output << "constrained_vars_ = " << it->first +1 << ";" << endl;
- output << "constrained_perfect_foresight_ = " << it->second.second << ";" << endl;
}
else
{
output << "constrained_vars_ = [constrained_vars_; " << it->first +1 << "];" << endl;
- output << "constrained_perfect_foresight_ = [constrained_perfect_foresight_; " << it->second.second << "];" << endl;
}
- const vector<AbstractShocksStatement::DetShockElement> &elems = it->second.first;
+ const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
for (int i = 0; i < (int) elems.size(); i++)
for (int j = elems[i].period1; j <= elems[i].period2; j++)
{
diff --git a/preprocessor/Shocks.hh b/preprocessor/Shocks.hh
index 539a4033634781dd059a10cdc967476f0581346e..68d8eee20a9ec83517c1e1d7deb66df118cb8ac4 100644
--- a/preprocessor/Shocks.hh
+++ b/preprocessor/Shocks.hh
@@ -41,7 +41,7 @@ public:
};
//The boolean element indicates if the shock is a surprise (false) or a perfect foresight (true) shock.
//This boolean is used only in case of conditional forecast with extended path method (simulation_type = deterministic).
- typedef map<int, pair< vector<DetShockElement>, bool> > det_shocks_t;
+ typedef map<int, vector<DetShockElement> > det_shocks_t;
protected:
//! Is this statement a "mshocks" statement ? (instead of a "shocks" statement)
const bool mshocks;