stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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).
unit-tests
Marco Ratto 2022-10-19 19:07:45 +02:00
parent 2ccd3d4a0e
commit d0150997f6
3 changed files with 86 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
matlab/+occbin/solve_one_constraint.m
View File

@ -125,20 +125,24 @@ for shock_period = 1:n_shocks_periods
regime_change_this_iteration=true; regime_change_this_iteration=true;
iter = 0; iter = 0;
nperiods_endogenously_increased = false;
guess_history_it = false; guess_history_it = false;
if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history
guess_history_it = true; guess_history_it = true;
end end
is_periodic=false; is_periodic=false;
is_periodic_loop=false;
binding_indicator_history={}; binding_indicator_history={};
max_err = NaN(max_iter,1); 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; iter = iter +1;
if binding_indicator(end) && nperiods_0<opts_simul_.max_periods if binding_indicator(end) && nperiods_0<opts_simul_.max_periods
binding_indicator = [binding_indicator; false ]; binding_indicator = [binding_indicator; false ];
nperiods_0 = nperiods_0 + 1; nperiods_0 = nperiods_0 + 1;
nperiods_endogenously_increased = true;
disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug) disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug)
elseif nperiods_0>=opts_simul_.max_periods elseif nperiods_0>=opts_simul_.max_periods
% enforce endogenously increased nperiods to not violate 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)); 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])); max_err(iter) = max(abs([err_viol;err_relax]));
regime_change_this_iteration = true; regime_change_this_iteration = true;
regime_violates_constraint_in_expectation(iter) = any(binding.constraint_1(1:end_periods) & ~binding_indicator(1:end_periods));
else else
regime_change_this_iteration = false; regime_change_this_iteration = false;
max_err(iter) = 0; 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); binding_indicator= (binding_indicator | binding.constraint_1) & ~(binding_indicator & relax.constraint_1);
end 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 for kiter=1:iter-1
vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)]; if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
is_periodic(kiter) = isequal(vvv, binding_indicator); % 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 end
is_periodic_all =is_periodic; is_periodic_all =is_periodic;
is_periodic = any(is_periodic); is_periodic = any(is_periodic);
if is_periodic && periodic_solution if is_periodic && periodic_solution
[merr,imerr]=min(max_err(find(is_periodic_all,1):end));
inx = find(is_periodic_all,1):iter; 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); inx = inx(imerr);
binding_indicator=binding_indicator_history{inx}; binding_indicator=binding_indicator_history{inx};
if inx<iter if inx<iter
@ -281,9 +313,14 @@ for shock_period = 1:n_shocks_periods
return return
end end
else 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) disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
if opts_simul_.waitbar; dyn_waitbar_close(hh); end
error_flag = 311; error_flag = 311;
end
if opts_simul_.waitbar; dyn_waitbar_close(hh); end
return return
end end
else else

47
matlab/+occbin/solve_two_constraints.m
View File

@ -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)); dyn_waitbar(shock_period/n_shocks_periods, hh, sprintf('Period %u of %u', shock_period,n_shocks_periods));
end end
regime_change_this_iteration=true; regime_change_this_iteration=true;
nperiods_endogenously_increased = false;
iter = 0; iter = 0;
guess_history_it = false; guess_history_it = false;
if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history if guess_history && (shock_period<=length(regime_history_guess)) %beyond guess regime history
guess_history_it = true; guess_history_it = true;
end end
is_periodic=false; is_periodic=false;
is_periodic_loop=false;
binding_indicator_history={}; binding_indicator_history={};
max_err = NaN(max_iter,1); 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; iter = iter +1;
if any(binding_indicator(end,:)) && nperiods_0<opts_simul_.max_periods if any(binding_indicator(end,:)) && nperiods_0<opts_simul_.max_periods
binding_indicator = [binding_indicator; false(1,2)]; binding_indicator = [binding_indicator; false(1,2)];
nperiods_0 = nperiods_0 + 1; nperiods_0 = nperiods_0 + 1;
nperiods_endogenously_increased = true;
disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug) disp_verbose(['nperiods has been endogenously increased up to ' int2str(nperiods_0) '.'],opts_simul_.debug)
elseif nperiods_0>=opts_simul_.max_periods elseif nperiods_0>=opts_simul_.max_periods
% enforce endogenously increased nperiods to not violate max_periods % enforce endogenously increased nperiods to not violate max_periods
@ -248,17 +252,43 @@ for shock_period = 1:n_shocks_periods
end end
binding_indicator = reshape(binding_indicator,nperiods_0+1,2); binding_indicator = reshape(binding_indicator,nperiods_0+1,2);
if iter>1 && regime_change_this_iteration if iter>1 && regime_change_this_iteration && ~nperiods_endogenously_increased
is_periodic=false(1,iter-1); % 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 for kiter=1:iter-1
vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 2)]; if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
is_periodic(kiter) = isequal(vvv, binding_indicator); % 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 end
is_periodic_all = is_periodic; is_periodic_all = is_periodic;
is_periodic = any(is_periodic); is_periodic = any(is_periodic);
if is_periodic && periodic_solution 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; 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); inx = inx(index_min_err);
binding_indicator=binding_indicator_history{inx}; %select regime history with same result, but smallest error binding_indicator=binding_indicator_history{inx}; %select regime history with same result, but smallest error
if inx<iter %update if needed if inx<iter %update if needed
@ -312,8 +342,13 @@ for shock_period = 1:n_shocks_periods
return; return;
end end
else 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) disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
error_flag = 311; error_flag = 311;
end
if opts_simul_.waitbar; dyn_waitbar_close(hh); end if opts_simul_.waitbar; dyn_waitbar_close(hh); end
return; return;
end end

2
matlab/get_error_message.m
View File

@ -188,6 +188,8 @@ switch info(1)
message = 'Occbin: Simulation did not converge, increase maxit or check_ahead_periods.'; message = 'Occbin: Simulation did not converge, increase maxit or check_ahead_periods.';
case 312 case 312
message = 'Occbin: Constraint(s) are binding at the end of the sample.'; 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 case 320
message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.'; message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.';
case 401 case 401