From d0150997f6f3b9ee909eaafedeaf27fefca0a36a Mon Sep 17 00:00:00 2001 From: Marco Ratto 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). --- 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 27d01a720..5b710cff5 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=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=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