diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index 7891373a183d4899fafb99edfc63504326984cd4..f7ae12636eabfdab6d17a2032616686bad5e2b09 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -70,7 +70,7 @@ while (t <= samplesize)
[endogenous_variables_paths(:,t), info_convergence, endogenousvariablespaths] = extended_path_core(ep.periods, DynareModel.endo_nbr, DynareModel.exo_nbr, innovations.positive_var_indx, ...
spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
DynareResults.steady_state, ...
- ep.verbosity, ep.use_bytecode, ep.stochastic.order, ...
+ verbosity, ep.use_bytecode, ep.stochastic.order, ...
DynareModel, pfm, ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
DynareOptions.lmmcp, ...
DynareOptions, ...
diff --git a/matlab/ep/extended_path_core.m b/matlab/ep/extended_path_core.m
index bc39e81fc06913e8ef62ee9d2687d6c771b7de7f..303e20660ce350e4214f847ec9cc5e6f5d6e53cb 100644
--- a/matlab/ep/extended_path_core.m
+++ b/matlab/ep/extended_path_core.m
@@ -59,42 +59,31 @@ if flag
options.stack_solve_algo = stack_solve_algo;
[tmp, maxerror] = perfect_foresight_solver_core(M, options, oo);
if maxerror>options.dynatol.f
- flag = false;
+ info_convergence = false;
else
- flag = true;
+ info_convergence = true;
end
- if ~flag && ~options.no_homotopy
- exo_orig = oo.exo_simul;
- endo_simul = repmat(steady_state,1,periods+1);
- for i = 1:10
- weight = i/10;
- oo.endo_simul = [weight*initial_conditions + (1-weight)*steady_state endo_simul];
- oo.exo_simul = repmat((1-weight)*oo.exo_steady_state', size(oo.exo_simul,1),1) + weight*exo_orig;
- [tmp, flag] = perfect_foresight_solver_core(M, options, oo);
- disp([i,flag])
- if ~flag
- break
- end
- endo_simul = tmp.endo_simul;
- end
- end
- info_convergence = flag;
else
switch(algo)
case 0
- [flag, endo_simul] = ...
+ [flag, tmp.endo_simul] = ...
solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1
- [flag, endo_simul] = ...
+ [flag, tmp.endo_simul] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
end
- tmp.endo_simul = endo_simul;
info_convergence = ~flag;
end
end
+if ~info_convergence && ~options.no_homotopy
+ [info_convergence, tmp.endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M, options, oo, pfm, ep, order, algo, 2, debug);
+end
+
if info_convergence
y = tmp.endo_simul(:,2);
else
y = NaN(size(endo_nbr,1));
end
+
+endogenousvariablespaths = tmp.endo_simul;
\ No newline at end of file
diff --git a/matlab/ep/extended_path_homotopy.m b/matlab/ep/extended_path_homotopy.m
new file mode 100644
index 0000000000000000000000000000000000000000..577f51ced202a7924d4e3ccafb7e390cae5846c7
--- /dev/null
+++ b/matlab/ep/extended_path_homotopy.m
@@ -0,0 +1,142 @@
+function [info_convergence, endo_simul] = extended_path_homotopy(endo_simul, exo_simul, M, options, oo, pfm, ep, order, algo, method, debug)
+
+% Copyright (C) 2016 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 ismember(method, [1, 2])
+ if isequal(method, 2)
+ endo_simul0 = endo_simul;
+ end
+ noconvergence = true;
+ iteration = 0;
+ weight = .1;
+ maxiter = 100;
+ increase_flag = false;
+ increase_factor = 1.2;
+ decrease_factor = 1.1;
+ state = false(5,1);
+ oldweight = weight;
+ while noconvergence
+ iteration = iteration + 1;
+ oo.endo_simul = endo_simul;
+ oo.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo.exo_steady_state));
+ if order==0
+ [tmp, maxerror] = perfect_foresight_solver_core(M, options, oo);
+ else
+ switch(algo)
+ case 0
+ [flag, tmp.endo_simul] = ...
+ solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
+ case 1
+ [flag, tmp.endo_simul] = ...
+ solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
+ end
+ end
+ if isequal(order, 0)
+ % Logical variable flag is false iff the solver fails.
+ flag = ~(maxerror>options.dynatol.f);
+ else
+ % Fix convention issue on the value of flag.
+ flag = ~flag;
+ end
+ if debug
+ disp(sprintf('%s\t %1.8f\t %s',int2str(iteration),weight,int2str(flag)))
+ end
+ if flag
+ state(2:end) = state(1:end-1);
+ state(1) = flag;
+ if isequal(weight, 1)
+ noconvergence = false;
+ break
+ end
+ if all(state)
+ increase_factor = 1+(increase_factor-1)*1.1;
+ state = false(size(state));
+ end
+ oldweight = weight;
+ weight = min(weight*increase_factor, 1);
+ increase_flag = true;
+ endo_simul = tmp.endo_simul;
+ else
+ state(2:end) = state(1:end-1);
+ state(1) = flag;
+ if increase_flag
+ weight = oldweight + (weight-oldweight)/100;
+ else
+ weight = min(weight/decrease_factor, 1);
+ end
+ end
+ if iteration>maxiter
+ break
+ end
+ if weight<1e-9
+ break
+ end
+ end
+ info_convergence = ~noconvergence;
+end
+
+if isequal(method, 3) || (isequal(method, 2) && noconvergence)
+ if isequal(method, 2)
+ endo_simul = endo_simul0;
+ end
+ weights = 0:(1/1000):1;
+ noconvergence = true;
+ index = 1;
+ jndex = 0;
+ nweights = length(weights);
+ while noconvergence
+ weight = weights(index);
+ oo.endo_simul = endo_simul;
+ oo.exo_simul = bsxfun(@plus, weight*exo_simul, (1-weight)*transpose(oo.exo_steady_state));
+ if order==0
+ [tmp, maxerror] = perfect_foresight_solver_core(M, options, oo);
+ else
+ switch(algo)
+ case 0
+ [flag, tmp.endo_simul] = ...
+ solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
+ case 1
+ [flag, tmp.endo_simul] = ...
+ solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
+ end
+ end
+ if isequal(order, 0)
+ % Logical variable flag is false iff the solver fails.
+ flag = ~(maxerror>options.dynatol.f);
+ else
+ % Fix convention issue on the value of flag.
+ flag = ~flag;
+ end
+ if debug
+ disp(sprintf('%s\t %1.8f\t %s',int2str(index),weight,int2str(flag)))
+ end
+ if flag
+ jndex = index;
+ if isequal(weight, 1)
+ noconvergence = false;
+ continue
+ end
+ index = index+1;
+ endo_simul = tmp.endo_simul;
+ else
+ break
+ end
+ end
+ info_convergence = ~noconvergence;
+end
\ No newline at end of file
diff --git a/matlab/ep/homotopic_steps.m b/matlab/ep/homotopic_steps.m
deleted file mode 100644
index c4839c46c052cbf5d9714f611c489f98371dc60e..0000000000000000000000000000000000000000
--- a/matlab/ep/homotopic_steps.m
+++ /dev/null
@@ -1,181 +0,0 @@
-function [info,tmp] = homotopic_steps(endo_simul0,exo_simul0,initial_weight,step_length,pfm)
-
-% Copyright (C) 2012 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_
-
-%Set bytecode flag
-bytecode_flag = options_.ep.use_bytecode;
-
-% Set increase and decrease factors.
-increase_factor = 5.0;
-decrease_factor = 0.2;
-
-% Save current state of oo_.endo_simul and oo_.exo_simul.
-endo_simul = endo_simul0;
-exxo_simul = exo_simul0;
-
-initial_step_length = step_length;
-max_iter = 1000/step_length;
-weight = initial_weight;
-verbose = options_.ep.verbosity;
-
-reduce_step_flag = 0;
-
-if verbose
- format long
- disp('Entering homotopic_steps')
-end
-
-% (re)Set iter.
-iter = 0;
-% (re)Set iter.
-jter = 0;
-% (re)Set weight.
-weight = initial_weight;
-% (re)Set exo_simul to zero.
-exo_simul0 = zeros(size(exo_simul0));
-while weight<1
- iter = iter+1;
- exo_simul0(2,:) = weight*exxo_simul(2,:);
- if bytecode_flag
- oo_.endo_simul = endo_simul_1;
- oo_.exo_simul = exo_simul_1;
- [flag,tmp] = bytecode('dynamic');
- else
- flag = 1;
- end
- if options_.debug
- save ep-test3 weight exo_simul0 endo_simul0
- end
- if flag
- if ~options_.ep.stochastic.order
- [flag,tmp,err] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
- else
- switch options_.ep.stochastic.algo
- case 0
- [flag,tmp] = ...
- solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
- case 1
- [flag,tmp] = ...
- solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_,pfm,options_.ep.stochastic.order,weight);
- % solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order);
- end
- end
- end
- info.convergence = ~flag;% Equal to one if the perfect foresight solver converged for the current value of weight.
- if verbose
- if info.convergence
- disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Ok!' ])
- else
- disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Convergence problem!' ])
- end
- end
- if info.convergence
- if options_.debug
- save ep-test2 weight exo_simul0 endo_simul0
- end
- endo_simul0 = tmp;
- jter = jter + 1;
- if jter>3
- if verbose
- disp('I am increasing the step length!')
- end
- step_length=step_length*increase_factor;
- jter = 0;
- end
- if abs(1-weight)<options_.dynatol.x;
- break
- end
- weight = weight+step_length;
- else% Perfect foresight solver failed for the current value of weight.
- if initial_weight>0 && abs(weight-initial_weight)<1e-12% First iteration, the initial weight is too high.
- if verbose
- disp('I am reducing the initial weight!')
- end
- initial_weight = initial_weight/2;
- weight = initial_weight;
- if weight<1e-12
- endo_simul0 = endo_simul;
- exo_simul0 = exxo_simul;
- info.convergence = 0;
- tmp = [];
- return
- end
- continue
- else% Initial weight is OK, but the perfect foresight solver failed on some subsequent iteration.
- if verbose
- disp('I am reducing the step length!')
- end
- jter = 0;
- if weight>0
- weight = weight-step_length;
- end
- step_length=step_length*decrease_factor;
- weight = weight+step_length;
- if step_length<options_.dynatol.x
- break
- end
- continue
- end
- end
- if iter>max_iter
- info = NaN;
- return
- end
-end
-if weight<1
- exo_simul0 = exxo_simul;
- if bytecode_flag
- oo_.endo_simul = endo_simul_1;
- oo_.exo_simul = exo_simul_1;
- [flag,tmp] = bytecode('dynamic');
- else
- flag = 1;
- end
- if flag
- if ~options_.ep.stochastic.order
- [flag,tmp,err] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
- else
- switch options_.ep.stochastic.algo
- case 0
- [flag,tmp] = ...
- solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
- case 1
- [flag,tmp] = ...
- solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_,pfm,options_.ep.stochastic.order,weight);
- % solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order);
- end
- end
- end
- info.convergence = ~flag;
- if info.convergence
- endo_simul0 = tmp;
- return
- else
- if step_length>options_.dynatol.x
- endo_simul0 = endo_simul;
- exo_simul0 = exxo_simul;
- info.convergence = 0;
- tmp = [];
- return
- else
- error('extended_path::homotopy: Oups! I did my best, but I am not able to simulate this model...')
- end
- end
-end