From fe1ba750eb0b61bfc1ff54c4184eabdc01b289b3 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Wed, 19 Oct 2022 19:07:45 +0200
Subject: [PATCH] accept periodic solution in simulations ONLY IF two regimes
differ by one period, to avoid pathological solutions.
We also do not check for periodicity when check ahead periods have been increased endogenously, again to avoid mis-identified periodicity.
Any other type of periodicity, is flagged as non-convergence with error code 313 (infinite loop of solutions).
(cherry picked from commit d0150997f6f3b9ee909eaafedeaf27fefca0a36a)
---
matlab/+occbin/solve_one_constraint.m | 49 ++++++++++++++++++++++----
matlab/+occbin/solve_two_constraints.m | 47 ++++++++++++++++++++----
matlab/get_error_message.m | 2 ++
3 files changed, 86 insertions(+), 12 deletions(-)
diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index 27d01a7204..5b710cff5a 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -125,20 +125,24 @@ for shock_period = 1:n_shocks_periods
regime_change_this_iteration=true;
iter = 0;
+ nperiods_endogenously_increased = false;
guess_history_it = false;
if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history
guess_history_it = true;
end
is_periodic=false;
+ is_periodic_loop=false;
binding_indicator_history={};
max_err = NaN(max_iter,1);
+ regime_violates_constraint_in_expectation = false(max_iter,1);
- while (regime_change_this_iteration && iter<max_iter && ~is_periodic)
+ while (regime_change_this_iteration && iter<max_iter && ~is_periodic && ~is_periodic_loop)
iter = iter +1;
if binding_indicator(end) && nperiods_0<opts_simul_.max_periods
binding_indicator = [binding_indicator; false ];
nperiods_0 = nperiods_0 + 1;
+ nperiods_endogenously_increased = true;
disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug)
elseif nperiods_0>=opts_simul_.max_periods
% enforce endogenously increased nperiods to not violate max_periods
@@ -201,6 +205,7 @@ for shock_period = 1:n_shocks_periods
err_relax = err.relax_constraint_1(relax.constraint_1(1:end_periods) & binding_indicator(1:end_periods));
max_err(iter) = max(abs([err_viol;err_relax]));
regime_change_this_iteration = true;
+ regime_violates_constraint_in_expectation(iter) = any(binding.constraint_1(1:end_periods) & ~binding_indicator(1:end_periods));
else
regime_change_this_iteration = false;
max_err(iter) = 0;
@@ -220,16 +225,43 @@ for shock_period = 1:n_shocks_periods
binding_indicator= (binding_indicator | binding.constraint_1) & ~(binding_indicator & relax.constraint_1);
end
- if iter>1 && regime_change_this_iteration
+ if iter>1 && regime_change_this_iteration && ~nperiods_endogenously_increased
+ % check for periodic solution only if nperiods is not
+ % increased endogenously
+ % first check for infinite loop
+ is_periodic_loop = false(iter-1,1);
for kiter=1:iter-1
- vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
- is_periodic(kiter) = isequal(vvv, binding_indicator);
+ if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+ % vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+ % is_periodic(kiter) = isequal(vvv, binding_indicator);
+ is_periodic_loop(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator);
+ else
+ is_periodic_loop(kiter) = false;
+ end
+ end
+ is_periodic_loop_all =is_periodic_loop;
+ is_periodic_loop = any(is_periodic_loop);
+ % only accept periodicity where regimes differ by one
+ % period!
+ is_periodic = false(iter-1,1);
+ for kiter=iter-1
+ if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+ % vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+ % is_periodic(kiter) = isequal(vvv, binding_indicator);
+ is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}-binding_indicator))==1;
+ else
+ is_periodic(kiter)=false;
+ end
end
is_periodic_all =is_periodic;
is_periodic = any(is_periodic);
if is_periodic && periodic_solution
- [merr,imerr]=min(max_err(find(is_periodic_all,1):end));
inx = find(is_periodic_all,1):iter;
+ inx1 = inx(find(~regime_violates_constraint_in_expectation(inx)));
+ if not(isempty(inx1))
+ inx=inx1;
+ end
+ [merr,imerr]=min(max_err(inx));
inx = inx(imerr);
binding_indicator=binding_indicator_history{inx};
if inx<iter
@@ -281,9 +313,14 @@ for shock_period = 1:n_shocks_periods
return
end
else
+ if is_periodic_loop
+ disp_verbose('Did not converge -- infinite loop of regimes.',opts_simul_.debug)
+ error_flag = 313;
+ else
disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
- if opts_simul_.waitbar; dyn_waitbar_close(hh); end
error_flag = 311;
+ end
+ if opts_simul_.waitbar; dyn_waitbar_close(hh); end
return
end
else
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index 4e448b241f..15c8ae1c2a 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -133,20 +133,24 @@ for shock_period = 1:n_shocks_periods
dyn_waitbar(shock_period/n_shocks_periods, hh, sprintf('Period %u of %u', shock_period,n_shocks_periods));
end
regime_change_this_iteration=true;
+ nperiods_endogenously_increased = false;
iter = 0;
guess_history_it = false;
if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history
guess_history_it = true;
end
is_periodic=false;
+ is_periodic_loop=false;
binding_indicator_history={};
max_err = NaN(max_iter,1);
+ regime_violates_constraint_in_expectation = false(max_iter,1);
- while (regime_change_this_iteration && iter<max_iter && ~is_periodic)
+ while (regime_change_this_iteration && iter<max_iter && ~is_periodic && ~is_periodic_loop)
iter = iter +1;
if any(binding_indicator(end,:)) && nperiods_0<opts_simul_.max_periods
binding_indicator = [binding_indicator; false(1,2)];
nperiods_0 = nperiods_0 + 1;
+ nperiods_endogenously_increased = true;
disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug)
elseif nperiods_0>=opts_simul_.max_periods
% enforce endogenously increased nperiods to not violate max_periods
@@ -248,17 +252,43 @@ for shock_period = 1:n_shocks_periods
end
binding_indicator = reshape(binding_indicator,nperiods_0+1,2);
- if iter>1 && regime_change_this_iteration
- is_periodic=false(1,iter-1);
+ if iter>1 && regime_change_this_iteration && ~nperiods_endogenously_increased
+ % check for periodic solution only if nperiods is not
+ % increased endogenously
+ % first check for infinite loop
+ is_periodic_loop = false(iter-1,1);
for kiter=1:iter-1
- vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 2)];
- is_periodic(kiter) = isequal(vvv, binding_indicator);
+ if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+ % vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+ % is_periodic(kiter) = isequal(vvv, binding_indicator);
+ is_periodic_loop(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator);
+ else
+ is_periodic_loop(kiter) = false;
+ end
+ end
+% is_periodic_loop_all =is_periodic_loop;
+ is_periodic_loop = any(is_periodic_loop);
+ % only accept periodicity where regimes differ by one
+ % period!
+ is_periodic=false(1,iter-1);
+ for kiter=iter-1
+ if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
+ % vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
+ % is_periodic(kiter) = isequal(vvv, binding_indicator);
+ is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}(:,1)-binding_indicator(:,1)))<=1 && length(find(binding_indicator_history{iter}(:,2)-binding_indicator(:,2)))<=1;
+ else
+ is_periodic(kiter)=false;
+ end
end
is_periodic_all = is_periodic;
is_periodic = any(is_periodic);
if is_periodic && periodic_solution
- [min_err,index_min_err]=min(max_err(find(is_periodic_all,1):end));
inx = find(is_periodic_all,1):iter;
+ inx1 = inx(find(~regime_violates_constraint_in_expectation(inx)));
+ if not(isempty(inx1))
+ inx=inx1;
+ end
+ [min_err,index_min_err]=min(max_err(inx));
inx = inx(index_min_err);
binding_indicator=binding_indicator_history{inx}; %select regime history with same result, but smallest error
if inx<iter %update if needed
@@ -312,8 +342,13 @@ for shock_period = 1:n_shocks_periods
return;
end
else
+ if is_periodic_loop
+ disp_verbose('Did not converge -- infinite loop of regimes.',opts_simul_.debug)
+ error_flag = 313;
+ else
disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
error_flag = 311;
+ end
if opts_simul_.waitbar; dyn_waitbar_close(hh); end
return;
end
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 1013184d69..cb236cfd2b 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -188,6 +188,8 @@ switch info(1)
message = 'Occbin: Simulation did not converge, increase maxit or check_ahead_periods.';
case 312
message = 'Occbin: Constraint(s) are binding at the end of the sample.';
+ case 313
+ message = 'Occbin: Simulation did not converge -- infinite loop of regimes';
case 320
message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.';
otherwise
--
GitLab