diff --git a/doc/dynare.texi b/doc/dynare.texi
index be979137002abd7327cbd6709b4d3da174fea2ca..3321ccd65abab69367259ce777b55b5a8c9a8a97 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -2123,6 +2123,11 @@ 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}
 
 For stochastic simulations, the @code{shocks} block specifies the non
@@ -4782,6 +4787,13 @@ 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} 
+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 rplot command.
 
 Finally, it is possible to do forecasting with a Bayesian VAR using
 the @code{bvar_forecast} command.
@@ -4941,7 +4953,10 @@ structural shocks that are needed to match the restricted paths. This
 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. Option
+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
 which will be matched to generate the constrained path.
 
@@ -4968,6 +4983,15 @@ 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
 
 @examplehead
@@ -4994,6 +5018,32 @@ conditional_forecast(parameter_set = calibration, controlled_varexo = (e, u), re
 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_forsight;
+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 ;
@@ -5006,6 +5056,11 @@ 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
 
 @deffn Command plot_conditional_forecast [@var{VARIABLE_NAME}@dots{}];
diff --git a/matlab/det_cond_forecast.m b/matlab/det_cond_forecast.m
new file mode 100644
index 0000000000000000000000000000000000000000..4decb95e653d9463146b1ca6915feaa020204cae
--- /dev/null
+++ b/matlab/det_cond_forecast.m
@@ -0,0 +1,352 @@
+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 = 60;
+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(exo_names(i,:), options_cond_fcst.controlled_varexo(j,:))
+            controlled_varexo(j) = i;
+        end
+    end
+end
+
+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_.solve_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 = controlled_varexo
+                if constrained_perfect_foresight(j)
+                    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;
+                else
+                    saved = oo_.exo_simul(maximum_lag+1,j);
+                    oo_.exo_simul(maximum_lag+1,j) = oo_.exo_simul(maximum_lag+1,j) + eps1;
+                    simul();
+                    J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
+                    oo_.exo_simul(maximum_lag+1,j) = 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 = controlled_varexo
+                    if constrained_perfect_foresight(j)
+                        oo_.exo_simul(time_index_constraint,j) = (oo_.exo_simul(time_index_constraint,j) + 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) = oo_.exo_simul(maximum_lag + 1,j) + 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 = controlled_varexo
+            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) = oo_.exo_simul(time,j);
+                end
+            else
+                exo(maximum_lag + t,j) = oo_.exo_simul(maximum_lag + 1,j);
+            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;
diff --git a/matlab/det_forecast.m b/matlab/det_forecast.m
new file mode 100644
index 0000000000000000000000000000000000000000..9e24fdb808af5e48cde164af778ca18767eb803f
--- /dev/null
+++ b/matlab/det_forecast.m
@@ -0,0 +1,163 @@
+function det_cond_forecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight)
+% Computes forecasts using the schocks retrieved from a condition forecast 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 = 60;
+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(exo_names(i,:), options_cond_fcst.controlled_varexo(j,:))
+            controlled_varexo(j) = i;
+        end
+    end
+end
+
+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_.solve_maxit;
+
+% Check the solution using a unconditional forecast (soft tune)
+
+initial_conditions = oo_.steady_state;
+terminal_conditions = oo_.steady_state;
+exo = oo_.exo_simul;
+T = options_.periods + 2;
+endo_simul = zeros(ny, T);
+endo_simul(:,1) = initial_conditions;
+endo_simul(:,T) = initial_conditions;
+exo_simul = zeros(T, nx);
+exo_simul(1,:) = [oo_.exo_steady_state' oo_.exo_det_steady_state'];
+exo_simul(T,:) = [oo_.exo_steady_state'  oo_.exo_det_steady_state'];
+past_val = 0;
+
+if pf && ~surprise
+    make_ex_;
+    make_y_;
+    oo_.endo_simul(:,1) = initial_conditions;
+    oo_.endo_simul(:,options_.periods + 2) = terminal_conditions;
+    %oo_.exo_simul = repmat(oo_.exo_steady_state, options_.periods + 2, 1);
+    oo_.exo_simul = exo;
+    simul();
+    endo_simul = oo_.endo_simul;
+    exo_simul = oo_.exo_simul;
+else
+    for t=1:constrained_periods
+        make_ex_;
+        make_y_;
+        disp(['t=' int2str(t) ' constrained_periods=' int2str(constrained_periods)]);
+        oo_.endo_simul(:,1) = initial_conditions;
+        oo_.endo_simul(:,options_.periods + 2) = terminal_conditions;
+        time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods - t + 1;
+        if t <= constrained_periods
+            for j = controlled_varexo
+                if constrained_perfect_foresight(j)
+                    for time = time_index_constraint;
+                        oo_.exo_simul(time,j) = exo(past_val + time,j);
+                    end;
+                    oo_.exo_simul(time+1, j)= oo_.exo_steady_state(j);
+                else
+                    oo_.exo_simul(maximum_lag + 1,j) = exo(maximum_lag + t,j);
+                end;
+            end;
+        else
+            tmp = size(oo_.exo_simul,1);
+            oo_.exo_simul = repmat(oo_.exo_steady_state',tmp,1);
+        end;
+        past_val = past_val + length(time_index_constraint);
+        simul();
+        initial_conditions = oo_.endo_simul(:,2);
+        if t < constrained_periods
+            endo_simul(:,t+1) = initial_conditions;
+            exo_simul(t+1,:) = oo_.exo_simul(2,:);
+        else
+            endo_simul(:,t + 1:t + options_cond_fcst.periods + maximum_lead) = oo_.endo_simul(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods + maximum_lead);
+            exo_simul(t+1,:) = oo_.exo_simul(2,:);
+        end;
+    end;
+end;
+oo_.endo_simul = endo_simul;
+oo_.exo_simul = exo_simul;
diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m
index 3a60857aabcb46827b1bf56cccbea491814b25f3..801ead636e0ff62b7880274d9d75782451ad5d30 100644
--- a/matlab/imcforecast.m
+++ b/matlab/imcforecast.m
@@ -1,4 +1,4 @@
-function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
+function imcforecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight)
 % Computes conditional forecasts.
 %
 % INPUTS
@@ -26,7 +26,7 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
 %  [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) 2006-2011 Dynare Team
+% Copyright (C) 2006-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -45,6 +45,15 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
 
 global options_ oo_ M_ bayestopt_
 
+
+if ~isempty(options_cond_fcst.simulation_type)
+    if strcmp(options_cond_fcst.simulation_type, 'deterministic')
+        disp('deterministic condtional forecast');
+        det_cond_forecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight);
+        return;
+    end
+end
+    
 if ~isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set)
     options_cond_fcst.parameter_set = 'posterior_mode';
 end
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 0c0e6578b9db2b3d53d1a97d0e1ec9f0d38140c5..cc755332570290282abdc1f7a351fa11f5ed08b4 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -1275,7 +1275,7 @@ void
 ConditionalForecastStatement::writeOutput(ostream &output, const string &basename) const
 {
   options_list.writeOutput(output, "options_cond_fcst_");
-  output << "imcforecast(constrained_paths_, constrained_vars_, options_cond_fcst_);" << endl;
+  output << "imcforecast(constrained_paths_, constrained_vars_, options_cond_fcst_, constrained_perfect_foresight_);" << endl;
 }
 
 PlotConditionalForecastStatement::PlotConditionalForecastStatement(int periods_arg, const SymbolList &symbol_list_arg) :
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index de89c232c9451834dd2eb1d67f727f7b791c20c1..546b65cedfed93fab69b52bff652912d59178869 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -96,7 +96,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 DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
+%token DATAFILE FILE DETERMINISTIC DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH
 %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
 %token <string_val> FLOAT_NUMBER
@@ -118,13 +118,14 @@ class ParsingDriver;
 %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS
 %token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF
 %token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED
-%token PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERIODS PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE
+%token PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERFECT_FORESIGHT PERIODS PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE
 %token PRINT PRIOR_MC PRIOR_TRUNC PRIOR_MODE PRIOR_MEAN POSTERIOR_MODE POSTERIOR_MEAN POSTERIOR_MEDIAN PRUNING
 %token <string_val> QUOTED_STRING
 %token QZ_CRITERIUM 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 SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO SOLVER_PERIODS
-%token STDERR STEADY STOCH_SIMUL SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME
+%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SIMULATION_TYPE
+%token SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL STOCHASTIC 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
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE
@@ -794,6 +795,12 @@ 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 ')'
@@ -864,8 +871,15 @@ check_options : steady_options;
 
 model_info : MODEL_INFO ';'
              { driver.model_info(); }
+            | MODEL_INFO '(' model_info_options_list ')' ';'
+              { driver.model_info(); }
            ;
 
+model_info_options_list : model_info_options_list COMMA model_info_options
+                   | model_info_options
+                   ;
+model_info_options :
+
 simul : SIMUL ';'
         { driver.simul(); }
       | SIMUL '(' simul_options_list ')' ';'
@@ -2112,6 +2126,7 @@ 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 ';'
@@ -2130,6 +2145,8 @@ 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(); }
@@ -2424,6 +2441,11 @@ 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_shocks : SHOCKS EQUAL '(' list_of_symbol_lists ')' { driver.option_symbol_list("shocks"); };
 o_labels : LABELS EQUAL '(' symbol_list ')' { driver.option_symbol_list("labels"); };
 o_ms_drop : DROP EQUAL INT_NUMBER { driver.option_num("ms.drop", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 3c8b40a9bc78220ff2de0a66ac450e9b2a8049ec..897597b1027e331674b8faf30df1cbe590dce5b2 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -458,6 +458,9 @@ string eofbuff;
 <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; }
@@ -481,6 +484,8 @@ string eofbuff;
 <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 e4e36f909cc1e728807ef18c3426ec1cae653387..d4520975fbd1c650d733eaa8403d506f4ba13803 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -606,10 +606,12 @@ ParsingDriver::add_det_shock(string *var, bool conditional_forecast)
       v.push_back(dse);
     }
 
-  det_shocks[symb_id] = v;
+  det_shocks[symb_id].first = v;
+  det_shocks[symb_id].second = det_shocks_expectation_pf;
 
   det_shocks_periods.clear();
   det_shocks_values.clear();
+  det_shocks_expectation_pf = false;
   delete var;
 }
 
@@ -749,6 +751,12 @@ 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 30ae85a40887c7e271b0fdd4c29692aaa4be27de..7c637ddc1cb1f1e2a5c8c6ad0fc126bcdc56174a 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -131,6 +131,8 @@ 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
@@ -332,6 +334,8 @@ 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 3536053e186c502145b214b1ef687858ddf099e4..31fcdd173c997813b0f1eba4fa938327585bacab 100644
--- a/preprocessor/Shocks.cc
+++ b/preprocessor/Shocks.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -46,11 +46,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.size(); i++)
+      for (size_t i = 0; i < it->second.first.size(); i++)
         {
-          const int &period1 = it->second[i].period1;
-          const int &period2 = it->second[i].period2;
-          const expr_t value = it->second[i].value;
+          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;
 
           if (period1 == period2)
             {
@@ -262,7 +262,7 @@ ConditionalForecastPathsStatement::checkPass(ModFileStructure &mod_file_struct,
        it != paths.end(); it++)
     {
       int this_path_length = 0;
-      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
+      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second.first;
       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);
@@ -289,11 +289,18 @@ ConditionalForecastPathsStatement::writeOutput(ostream &output, const string &ba
        it != paths.end(); it++)
     {
       if (it == paths.begin())
-        output << "constrained_vars_ = " << it->first +1 << ";" << endl;
+        {
+          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_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;
+      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second.first;
       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 7eda7f7aedc6ac4e326407dba3f2282c09a22aed..de20945c72861c82e3764aa91cb7dd8fbadce2a6 100644
--- a/preprocessor/Shocks.hh
+++ b/preprocessor/Shocks.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -39,7 +39,9 @@ public:
     int period2;
     expr_t value;
   };
-  typedef map<int, vector<DetShockElement> > det_shocks_t;
+  //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;
 protected:
   //! Is this statement a "mshocks" statement ? (instead of a "shocks" statement)
   const bool mshocks;