diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index 056b5f4ed..91d4def6b 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -118,10 +118,12 @@ replic_nbr = ep.replic_nbr;
switch ep.innovation_distribution
case 'gaussian'
shocks = transpose(transpose(covariance_matrix_upper_cholesky)* ...
- randn(effective_number_of_shocks,sample_size*replic_nbr));
+ randn(effective_number_of_shocks,sample_size* ...
+ replic_nbr));
+ shocks(:,positive_var_indx) = shocks;
case 'calibrated'
replic_nbr = 1;
- shocks = zeros(sample_size,effective_number_of_shocks);
+ shocks = zeros(sample_size,M_.exo_nbr);
for i = 1:length(M_.unanticipated_det_shocks)
k = M_.unanticipated_det_shocks(i).periods;
ivar = M_.unanticipated_det_shocks(i).exo_id;
@@ -132,7 +134,6 @@ switch ep.innovation_distribution
socks(k,ivar) = shocks(k,ivar) * v;
end
end
- shocks = shocks(:,positive_var_indx);
otherwise
error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
@@ -203,7 +204,7 @@ while (t <= sample_size)
parfor k = 1:replic_nbr
exo_simul = repmat(oo_.exo_steady_state',periods+2,1);
% exo_simul(1:sample_size+3-t,:) = exo_simul_(t:end,:);
- exo_simul(2,positive_var_indx) = exo_simul_(M_.maximum_lag+t,positive_var_indx) + ...
+ exo_simul(2,:) = exo_simul_(M_.maximum_lag+t,:) + ...
shocks((t-2)*replic_nbr+k,:);
initial_conditions = results{k}(:,t-1);
results{k}(:,t) = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
@@ -219,8 +220,8 @@ while (t <= sample_size)
maximum_lead,1);
% exo_simul(1:sample_size+maximum_lag+maximum_lead-t+1,:) = ...
% exo_simul_(t:end,:);
- exo_simul(maximum_lag+1,positive_var_indx) = ...
- exo_simul_(maximum_lag+t,positive_var_indx) + shocks((t-2)*replic_nbr+k,:);
+ exo_simul(maximum_lag+1,:) = ...
+ exo_simul_(maximum_lag+t,:) + shocks((t-2)*replic_nbr+k,:);
initial_conditions = results{k}(:,t-1);
results{k}(:,t) = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
exo_simul,ep.init,initial_conditions,...
diff --git a/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m b/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
index 3df0a788b..286bd16ac 100644
--- a/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
+++ b/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
@@ -20,7 +20,9 @@ function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,Dynar
pfm.lead_lag_incidence = DynareModel.lead_lag_incidence;
pfm.ny = DynareModel.endo_nbr;
pfm.Sigma = DynareModel.Sigma_e;
-pfm.Omega = chol(pfm.Sigma,'upper'); % Sigma = Omega'*Omega
+if det(pfm.Sigma) > 0
+ pfm.Omega = chol(pfm.Sigma,'upper'); % Sigma = Omega'*Omega
+end
pfm.number_of_shocks = length(pfm.Sigma);
pfm.stochastic_order = DynareOptions.ep.stochastic.order;
pfm.max_lag = DynareModel.maximum_endo_lag;
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index b4d8a632a..7cbe7f654 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -93,8 +93,8 @@ o_periods = periods;
ZERO = zeros(length(i_upd),1);
h1 = clock ;
-iA = zeros(periods*M_.NNZDerivatives(1),3);
-for iter = 1:options_.simul.maxit
+iA = zeros(periods*M.NNZDerivatives(1),3);
+for iter = 1:options.simul.maxit
h2 = clock ;
i_rows = (1:ny)';
diff --git a/matlab/perfect_foresight_solver.m b/matlab/perfect_foresight_solver.m
deleted file mode 100644
index 0e03d9646..000000000
--- a/matlab/perfect_foresight_solver.m
+++ /dev/null
@@ -1,160 +0,0 @@
-function perfect_foresight_solver()
-% Computes deterministic simulations
-%
-% INPUTS
-% None
-%
-% OUTPUTS
-% none
-%
-% ALGORITHM
-%
-% SPECIAL REQUIREMENTS
-% none
-
-% Copyright (C) 1996-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 .
-
-global M_ options_ oo_
-
-if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 7
- error('PERFECT_FORESIGHT_SOLVER: stack_solve_algo must be between 0 and 7')
-end
-
-if ~options_.block && ~options_.bytecode && options_.stack_solve_algo ~= 0 ...
- && options_.stack_solve_algo ~= 6 && options_.stack_solve_algo ~= 7
- error('PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo=0 or stack_solve_algo=6 when not using block nor bytecode option')
-end
-
-if options_.block && ~options_.bytecode && options_.stack_solve_algo == 5
- error('PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 5 without bytecode option')
-end
-
-if (options_.block || options_.bytecode) && options_.stack_solve_algo == 6
- error('PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 6 with block or bytecode option')
-end
-
-if isoctave && options_.stack_solve_algo == 2
- error('PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 2 under Octave')
-end
-
-
-if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maximum_lag+options_.periods+M_.maximum_lead ])
- error('PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
-end
-
-if isempty(oo_.exo_simul) || any(size(oo_.exo_simul) ~= [ M_.maximum_lag+options_.periods+M_.maximum_lead, M_.exo_nbr ])
- error('PERFECT_FORESIGHT_SOLVER: ''oo_.exo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
-end
-
-
-if isempty(options_.scalv) || options_.scalv == 0
- options_.scalv = oo_.steady_state;
-end
-
-options_.scalv= 1;
-
-if options_.debug
- model_static = str2func([M_.fname,'_static']);
- for ii=1:size(oo_.exo_simul,1)
- [residual(:,ii)] = model_static(oo_.steady_state, oo_.exo_simul(ii,:),M_.params);
- end
- problematic_periods=find(any(isinf(residual)) | any(isnan(residual)))-M_.maximum_endo_lag;
- if ~isempty(problematic_periods)
- period_string=num2str(problematic_periods(1));
- for ii=2:length(problematic_periods)
- period_string=[period_string, ', ', num2str(problematic_periods(ii))];
- end
- fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string);
- fprintf('WARNING: Check for division by 0.\n')
- end
-end
-
-% Effectively compute simulation, possibly with homotopy
-if options_.no_homotopy
- [oo_.endo_simul,oo_.deterministic_simulation.status] = perfect_foresight_solver_core(M_,oo_,options_);
-else
- exosim = oo_.exo_simul;
- exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
- endosim = oo_.endo_simul;
- endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead);
-
- current_weight = 0; % Current weight of target point in convex combination
- step = 1;
- success_counter = 0;
-
- while (step > options_.dynatol.x)
-
- new_weight = current_weight + step; % Try this weight, and see if it succeeds
- if new_weight >= 1
- new_weight = 1; % Don't go beyond target point
- step = new_weight - current_weight;
- end
-
- % Compute convex combination for exo path and initial/terminal endo conditions
- % But take care of not overwriting the computed part of oo_.endo_simul
- oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight);
- endocombi = endosim*new_weight + endoinit*(1-new_weight);
- oo_.endo_simul(:,1:M_.maximum_endo_lag) = endocombi(:,1:M_.maximum_endo_lag);
- oo_.endo_simul(:,(end-M_.maximum_endo_lead):end) = endocombi(:,(end-M_.maximum_endo_lead):end);
-
- saved_endo_simul = oo_.endo_simul;
-
- [oo_.endo_simul,oo_.deterministic_simulation.status] = perfect_foresight_solver_core(M_,oo_,options_);
-
- if oo_.deterministic_simulation.status == 1
- current_weight = new_weight;
- if current_weight >= 1
- break
- end
- success_counter = success_counter + 1;
- if success_counter >= 3
- success_counter = 0;
- step = step * 2;
- disp([ 'Homotopy step succeeded, doubling step size (completed ' sprintf('%.1f', current_weight*100) '%, step size ' sprintf('%.3g', step) ')' ])
- else
- disp([ 'Homotopy step succeeded (completed ' sprintf('%.1f', current_weight*100) '%, step size ' sprintf('%.3g', step) ')' ])
- end
- else
- oo_.endo_simul = saved_endo_simul;
- success_counter = 0;
- step = step / 2;
- disp([ 'Homotopy step failed, halving step size (completed ' sprintf('%.1f', current_weight*100) '%, step size ' sprintf('%.3g', step) ')' ])
- end
- end
-end
-
-if oo_.deterministic_simulation.status == 1
- disp('Perfect foresight solution found.')
-else
- warning('Failed to solve perfect foresight model')
-end
-
-dyn2vec;
-
-if isnan(options_.initial_period)
- initial_period = dates(1,1);
-else
- initial_period = options_.initial_period;
-end
-
-ts = dseries(transpose(oo_.endo_simul),initial_period,cellstr(M_.endo_names));
-assignin('base', 'Simulated_time_series', ts);
-
-end
-
-
diff --git a/matlab/perfect_foresight_solver_core.m b/matlab/perfect_foresight_solver_core.m
deleted file mode 100644
index 3b2119130..000000000
--- a/matlab/perfect_foresight_solver_core.m
+++ /dev/null
@@ -1,134 +0,0 @@
-function [endo_simul, status] = perfect_foresight_solver_core(M,oo,options)
-% Core function to compute deterministic simulations
-%
-% INPUTS
-% M: model structure
-% oo: output structure
-% options: options structure
-%
-% OUTPUTS
-% endo_simul: matrix endogenous variables
-% deterministic_simulation: simulation status
-%
-% ALGORITHM
-%
-% various
-%
-% SPECIAL REQUIREMENTS
-% none
-
-% Copyright (C) 1996-2015 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 .
-
-endo_simul = oo.endo_simul;
-status = 0;
-deterministic_simulation = struct();
-
-if(options.block)
- if(options.bytecode)
- [info, endo_simul] = bytecode('dynamic');
- if info == 1
- status = 0;
- else
- status = 1;
- end
- mexErrCheck('bytecode', info);
- else
- eval([M.fname '_dynamic']);
- end
-else
- if(options.bytecode)
- [info, endo_simul]=bytecode('dynamic');
- if info == 1
- status = 0;
- else
- status = 1;
- end;
- mexErrCheck('bytecode', info);
- else
- if M.maximum_endo_lead == 0
- % Purely backward model
- global oo_ options_
- oo_ = oo;
- options_ = options;
- sim1_purely_backward;
- endo_simul = oo_.endo_simul;
- if oo_.deterministic_simulation.status == 1
- status = 0;
- end
- elseif M.maximum_endo_lag == 0
- % Purely forward model
- global oo_ options_
- oo_ = oo;
- options_ = options;
- sim1_purely_forward;
- endo_simul = oo_.endo_simul;
- if oo_.deterministic_simulation.status == 1
- status = 0;
- end
- else
- % General case
- if options.stack_solve_algo == 0
- oo = sim1(M,options,oo);
- endo_simul = oo.endo_simul;
- if oo.deterministic_simulation.status == 1
- status = 0;
- end
- elseif options.stack_solve_algo == 6
- global oo_ options_
- oo_ = oo;
- options_ = options;
- sim1_lbj;
- endo_simul = oo_.endo_simul;
- if oo_.deterministic_simulation.status == 1
- status = 0;
- end
- elseif options.stack_solve_algo == 7
- periods = options.periods;
- if ~isfield(options.lmmcp,'lb')
- [lb,ub,pfm.eq_index] = get_complementarity_conditions(M);
- options.lmmcp.lb = repmat(lb,periods,1);
- options.lmmcp.ub = repmat(ub,periods,1);
- end
-
- y = endo_simul;
- y0 = y(:,1);
- yT = y(:,periods+2);
- z = y(:,2:periods+1);
- illi = M.lead_lag_incidence';
- [i_cols,~,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,~,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M.lead_lag_incidence(1:2,:)');
- [y,info] = dynare_solve(@perfect_foresight_problem,z(:),1, ...
- str2func([M.fname '_dynamic']),y0,yT, ...
- oo.exo_simul,M.params,oo.steady_state, ...
- options.periods,M.endo_nbr,i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
- M.NNZDerivatives(1));
- endo_simul = [y0 reshape(y,M.endo_nbr,periods) yT];
- if info == 1
- status = 1;
- else
- status = 0;
- end;
- end
- end
- end
-end
-
-