Skip to content
Snippets Groups Projects
Commit 3fa6a96b authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Added homotopy. 

Fixes the issue of paths containing spurious solutions for (stochastic)
perfect foresight models.
parent cfc69576
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
matlab/ep/extended_path.m
+ 1
− 1
View file @ 3fa6a96b
...@@ -70,7 +70,7 @@ while (t <= samplesize) ...@@ -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, ... [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), ... spfm_exo_simul, ep.init, endogenous_variables_paths(:,t-1), ...
DynareResults.steady_state, ... 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, ... DynareModel, pfm, ep.stochastic.algo, ep.solve_algo, ep.stack_solve_algo, ...
DynareOptions.lmmcp, ... DynareOptions.lmmcp, ...
DynareOptions, ... DynareOptions, ...
......
matlab/ep/extended_path_core.m
+ 10
− 21
View file @ 3fa6a96b
...@@ -59,42 +59,31 @@ if flag ...@@ -59,42 +59,31 @@ if flag
options.stack_solve_algo = stack_solve_algo; options.stack_solve_algo = stack_solve_algo;
[tmp, maxerror] = perfect_foresight_solver_core(M, options, oo); [tmp, maxerror] = perfect_foresight_solver_core(M, options, oo);
if maxerror>options.dynatol.f if maxerror>options.dynatol.f
flag = false; info_convergence = false;
else 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 end
info_convergence = flag;
else else
switch(algo) switch(algo)
case 0 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); solve_stochastic_perfect_foresight_model(endo_simul, exo_simul, pfm, ep.stochastic.quadrature.nodes, ep.stochastic.order);
case 1 case 1
[flag, endo_simul] = ... [flag, tmp.endo_simul] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order); solve_stochastic_perfect_foresight_model_1(endo_simul, exo_simul, options, pfm, ep.stochastic.order);
end end
tmp.endo_simul = endo_simul;
info_convergence = ~flag; info_convergence = ~flag;
end end
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 if info_convergence
y = tmp.endo_simul(:,2); y = tmp.endo_simul(:,2);
else else
y = NaN(size(endo_nbr,1)); y = NaN(size(endo_nbr,1));
end end
endogenousvariablespaths = tmp.endo_simul;
\ No newline at end of file
matlab/ep/extended_path_homotopy.m 0 → 100644
+ 142
− 0
View file @ 3fa6a96b
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
matlab/ep/homotopic_steps.m deleted 100644 → 0
+ 0
− 181
View file @ cfc69576
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment