Skip to content
GitLab
Commit 192b5c71 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Fixed bugs in routines related to the extended path algorithm.

parent 0aaa7a8a
Hide whitespace changes
Inline Side-by-side
matlab/ep/extended_path.m
View file @ 192b5c71
......@@ -105,6 +105,9 @@ periods = options_.ep.periods;
pfm.periods = options_.ep.periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
% keep a copy of pfm.i_upd
i_upd = pfm.i_upd;
% Set the algorithm for the perfect foresight solver
options_.stack_solve_algo = options_.ep.stack_solve_algo;
......@@ -190,23 +193,27 @@ while (t<sample_size)
initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
else
endo_simul_1 = repmat(steady_state,1,periods1+2);
if t==1
endo_simul_1 = repmat(steady_state,1,periods1+2);
end
end
% Solve a perfect foresight model.
increase_periods = 0;
% Keep a copy of endo_simul_1
endo_simul = endo_simul_1;
while 1
if ~increase_periods
if bytecode_flag
oo_.endo_simul = endo_simul_1;
oo_.exo_simul = exo_simul_1;
[flag,tmp] = bytecode('dynamic');
flag
pause
else
flag = 1;
end
if flag
[flag,tmp] = solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,...
pfm1, ...
options_.ep.nnodes,...
options_.ep.order);
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
end
info_convergence = ~flag;
end
......@@ -275,6 +282,8 @@ while (t<sample_size)
end
% Solve the perfect foresight model with an increased number of periods.
if bytecode_flag
oo_.endo_simul = endo_simul_1;
oo_.exo_simul = exo_simul_1;
[flag,tmp] = bytecode('dynamic');
else
flag = 1;
......@@ -327,14 +336,30 @@ while (t<sample_size)
end
end% while
if ~info_convergence% If exited from the while loop without achieving convergence, use an homotopic approach
[INFO,tmp] = homotopic_steps(.5,.01,pfm1);
if (~isstruct(INFO) && isnan(INFO)) || ~info_convergence
[INFO,tmp] = homotopic_steps(0,.01,pfm1);
if ~do_not_check_stability_flag
periods1 = ep.periods;
pfm1.periods = periods1;
pfm1.i_upd = i_upd;
exo_simul_1 = exo_simul_1(1:(periods1+2),:);
endo_simul_1 = endo_simul_1(:,1:(periods1+2));
end
[INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,.5,.01,pfm);
if isstruct(INFO)
info_convergence = INFO.convergence;
else
info_convergence = 0;
end
if ~info_convergence
[INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,0,.01,pfm);
if isstruct(INFO)
info_convergence = INFO.convergence;
else
info_convergence = 0;
end
if ~info_convergence
disp('Homotopy:: No convergence of the perfect foresight model solver!')
error('I am not able to simulate this model!');
else
info_convergence = 1;
endo_simul_1 = tmp;
if verbosity && info_convergence
disp('Homotopy:: Convergence of the perfect foresight model solver!')
......@@ -349,7 +374,10 @@ while (t<sample_size)
end
end
% Save results of the perfect foresight model solver.
time_series(:,t) = tmp;
time_series(:,t) = endo_simul_1(:,2);
endo_simul_1(:,1:end-1) = endo_simul_1(:,2:end);
endo_simul_1(:,1) = time_series(:,t);
endo_simul_1(:,end) = oo_.steady_state;
end% (while) loop over t
dyn_waitbar_close(hh);
......
matlab/ep/homotopic_steps.m
View file @ 192b5c71
function [info,tmp] = homotopic_steps(initial_weight,step_length,pfm)
global oo_ options_ M_
function [info,tmp] = homotopic_steps(endo_simul0,exo_simul0,initial_weight,step_length,pfm)
global options_ oo_
%Set bytecode flag
bytecode_flag = options_.ep.use_bytecode;
......@@ -9,26 +9,8 @@ increase_factor = 5.0;
decrease_factor = 0.2;
% Save current state of oo_.endo_simul and oo_.exo_simul.
endo_simul = oo_.endo_simul;
exxo_simul = oo_.exo_simul;
[idx,jdx] = find(abs(exxo_simul)>1e-12);
idx = unique(idx);
if idx(1)-2
% The first scalar in idx must be equal to two.
error('extended_path::homotopy: Something is wrong in oo_.exo_simul!')
end
stochastic_extended_path_depth = length(idx);
if stochastic_extended_path_depth>1
for i=2:length(idx)
if (idx(i)-idx(i-1)-1)
error('extended_path::homotopy: Something is wrong in oo_.exo_simul!')
end
end
end
endo_simul = endo_simul0;
exxo_simul = exo_simul0;
initial_step_length = step_length;
max_iter = 1000/step_length;
......@@ -41,115 +23,114 @@ if verbose
format long
end
for d=1:stochastic_extended_path_depth
% (re)Set iter.
iter = 0;
% (re)Set iter.
jter = 0;
% (re)Set weight.
weight = initial_weight;
% (re)Set exo_simul to zero.
oo_.exo_simul = zeros(size(oo_.exo_simul));
while weight<1
iter = iter+1;
oo_.exo_simul(idx(d),:) = weight*exxo_simul(idx(d),:);
if d>1
oo_.exo_simul(1:idx(d-1),:) = exxo_simul(1:idx(d-1),:);
end
if bytecode_flag
[flag,tmp] = bytecode('dynamic');
% (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 flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
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
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(oo_.endo_simul,oo_.exo_simul,pfm);
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Convergence problem!' ])
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
if info.convergence
%if d<stochastic_extended_path_depth
endo_simul0 = tmp;
%end
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 info.convergence
if d<stochastic_extended_path_depth
oo_.endo_simul = tmp;
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
jter = jter + 1;
if jter>3
if verbose
disp('I am increasing the step length!')
end
step_length=step_length*increase_factor;
jter = 0;
initial_weight = initial_weight/2;
weight = initial_weight;
if weight<1e-12
endo_simul0 = endo_simul;
exo_simul0 = exxo_simul;
info.convergence = 0;
info.depth = d;
tmp = [];
return
end
if abs(1-weight)<options_.dynatol.x;
break
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;
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
oo_.endo_simul = endo_simul;
oo_.exo_simul = exxo_simul;
info.convergence = 0;
info.depth = d;
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
if step_length<options_.dynatol.x
break
end
end
if iter>max_iter
info = NaN;
return
continue
end
end
if weight<1
oo_.exo_simul = exxo_simul;
if bytecode_flag
[flag,tmp] = bytecode('dynamic');
else
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(oo_.endo_simul,oo_.exo_simul,pfm);
end
info.convergence = ~flag;
if info.convergence
oo_.endo_simul = tmp;
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
[flag,tmp] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
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;
info.depth = d;
tmp = [];
return
else
if step_length>options_.dynatol.x
oo_.endo_simul = endo_simul;
oo_.exo_simul = exxo_simul;
info.convergence = 0;
info.depth = d;
tmp = [];
return
else
error('extended_path::homotopy: Oups! I did my best, but I am not able to simulate this model...')
end
error('extended_path::homotopy: Oups! I did my best, but I am not able to simulate this model...')
end
end
end
\ No newline at end of file
end
matlab/solve_perfect_foresight_model.m
View file @ 192b5c71
......@@ -3,6 +3,7 @@ function [flag,endo_simul,err] = solve_perfect_foresight_model(endo_simul,exo_si
flag = 0;
err = 0;
stop = 0;
nan_flag = 0;
model_dynamic = pfm.dynamic_model;
......@@ -57,6 +58,10 @@ function [flag,endo_simul,err] = solve_perfect_foresight_model(endo_simul,exo_si
break
end
dy = -A\res;
if any(isnan(dy))
nan_flag = 1;
break
end
Y(pfm.i_upd) = Y(pfm.i_upd) + dy;
end
......@@ -71,6 +76,17 @@ function [flag,endo_simul,err] = solve_perfect_foresight_model(endo_simul,exo_si
flag = 1;% more iterations are needed.
endo_simul = 1;
end
if nan_flag
if pfm.verbose
fprintf('\n') ;
disp([' Total time of simulation :' num2str(etime(clock,h1))]) ;
fprintf('\n') ;
disp(['WARNING : NaNs!']) ;
fprintf('\n') ;
end
flag = 1;
endo_simul = 1;
end
if pfm.verbose
disp (['-----------------------------------------------------']) ;
end
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment