diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m index 9e02d4a5f90199f0db4e28543803a02d77a93f45..1bcd697e23bd283db176fff1f288cc648f158cd7 100644 --- a/matlab/ep/extended_path.m +++ b/matlab/ep/extended_path.m @@ -248,94 +248,84 @@ while (t<sample_size) end end end - % Test if periods is big enough. The variable delta measures the maximum absolute variation during - % the last periods of the simulated path. This variation has to be close to zero (because the - % economy is assumed to be in the steady state at the end of the simulated path). - if info.convergence - delta = max(max(abs(tmp(:,end-options_.ep.lp:end)-tmp(:,end-options_.ep.lp-1:end-1)))); - end - if info.convergence && ~increase_periods && delta<options_.dynatol.x - % Exit from the while loop (the number of periods is big enough). - break - else - % Increase the number of periods. - options_.periods = options_.periods + options_.ep.step; - pfm.periods = options_.periods; - pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); - % Increment the counter. - increase_periods = increase_periods + 1; - if verbosity - if t<10 - disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) - elseif t<100 - disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) - elseif t<1000 - disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) - else - disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) - end - end - if info.convergence - % If the previous call to the perfect foresight model solver exited - % announcing that the routine converged, adapt the size of oo_.endo_simul - % and oo_.exo_simul. - oo_.endo_simul = [ tmp , repmat(oo_.steady_state,1,options_.ep.step) ]; - oo_.exo_simul = [ oo_.exo_simul ; zeros(options_.ep.step,size(shocks,2)) ]; - tmp_old = tmp; + % Test if periods is big enough. + % Increase the number of periods. + options_.periods = options_.periods + options_.ep.step; + pfm.periods = options_.periods; + pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); + % Increment the counter. + increase_periods = increase_periods + 1; + if verbosity + if t<10 + disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) + elseif t<100 + disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) + elseif t<1000 + disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) else - % If the previous call to the perfect foresight model solver exited - % announcing that the routine did not converge, then tmp=1... Maybe - % should change that, because in some circonstances it may usefull - % to know where the routine did stop, even if convergence was not - % achieved. - oo_.endo_simul = [ oo_.endo_simul , repmat(oo_.steady_state,1,options_.ep.step) ]; - oo_.exo_simul = [ oo_.exo_simul ; zeros(options_.ep.step,size(shocks,2)) ]; + disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(options_.periods) '.']) end - % Solve the perfect foresight model with an increased number of periods. - [flag,tmp] = bytecode('dynamic'); - if flag - [flag,tmp] = solve_perfect_foresight_model(oo_.endo_simul,oo_.exo_simul,pfm); + end + if info.convergence + % If the previous call to the perfect foresight model solver exited + % announcing that the routine converged, adapt the size of oo_.endo_simul + % and oo_.exo_simul. + oo_.endo_simul = [ tmp , repmat(oo_.steady_state,1,options_.ep.step) ]; + oo_.exo_simul = [ oo_.exo_simul ; zeros(options_.ep.step,size(shocks,2)) ]; + tmp_old = tmp; + else + % If the previous call to the perfect foresight model solver exited + % announcing that the routine did not converge, then tmp=1... Maybe + % should change that, because in some circonstances it may usefull + % to know where the routine did stop, even if convergence was not + % achieved. + oo_.endo_simul = [ oo_.endo_simul , repmat(oo_.steady_state,1,options_.ep.step) ]; + oo_.exo_simul = [ oo_.exo_simul ; zeros(options_.ep.step,size(shocks,2)) ]; + end + % Solve the perfect foresight model with an increased number of periods. + [flag,tmp] = bytecode('dynamic'); + if flag + [flag,tmp] = solve_perfect_foresight_model(oo_.endo_simul,oo_.exo_simul,pfm); + end + info.convergence = ~flag; + if info.convergence + % If the solver achieved convergence, check that simulated paths did not + % change during the first periods. + % Compute the maximum deviation between old path and new path over the + % first periods + delta = max(max(abs(tmp(:,2:options_.ep.fp)-tmp_old(:,2:options_.ep.fp)))); + if delta<options_.dynatol.x + % If the maximum deviation is close enough to zero, reset the number + % of periods to ep.periods + options_.periods = options_.ep.periods; + pfm.periods = options_.periods; + pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); + % Cut oo_.exo_simul and oo_.endo_simul consistently with the resetted + % number of periods and exit from the while loop. + oo_.exo_simul = oo_.exo_simul(1:(options_.periods+2),:); + oo_.endo_simul = oo_.endo_simul(:,1:(options_.periods+2)); + break end - info.convergence = ~flag; - if info.convergence - % If the solver achieved convergence, check that simulated paths did not - % change during the first periods. - % Compute the maximum deviation between old path and new path over the - % first periods - delta = max(max(abs(tmp(:,2:options_.ep.fp)-tmp_old(:,2:options_.ep.fp)))); - if delta<options_.dynatol.x - % If the maximum deviation is close enough to zero, reset the number - % of periods to ep.periods - options_.periods = options_.ep.periods; - pfm.periods = options_.periods; - pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); - % Cut oo_.exo_simul and oo_.endo_simul consistently with the resetted - % number of periods and exit from the while loop. - oo_.exo_simul = oo_.exo_simul(1:(options_.periods+2),:); - oo_.endo_simul = oo_.endo_simul(:,1:(options_.periods+2)); - break - end - else - % The solver did not converge... Try to solve the model again with a bigger - % number of periods, except if the number of periods has been increased more - % than 10 times. - if increase_periods==10; - if verbosity - if t<10 - disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) - elseif t<100 - disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) - elseif t<1000 - disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) - else - disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) - end + else + % The solver did not converge... Try to solve the model again with a bigger + % number of periods, except if the number of periods has been increased more + % than 10 times. + if increase_periods==10; + if verbosity + if t<10 + disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) + elseif t<100 + disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) + elseif t<1000 + disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) + else + disp(['Time: ' int2str(t) '. Even with ' int2str(options_.periods) ', I am not able to solve the perfect foresight model. Use homotopy instead...']) end - % Exit from the while loop. - break end - end% if info.convergence - end% ~increase_periods && delta<options_.dynatol.x + % Exit from the while loop. + break + end + end% if info.convergence end% while if ~info.convergence% If exited from the while loop without achieving convergence, use an homotopic approach [INFO,tmp] = homotopic_steps(.5,.01);