diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m
index 5bd5f9098dc2dcce839f09a32d3e6e9781aaf01d..384fa57bdf8be48c24c3f2dfdd9618b43b01ff6d 100644
--- a/matlab/+occbin/DSGE_smoother.m
+++ b/matlab/+occbin/DSGE_smoother.m
@@ -106,6 +106,7 @@ opts_simul.maxit = options_.occbin.smoother.maxit;
opts_simul.waitbar = options_.occbin.smoother.waitbar;
opts_simul.periods = options_.occbin.smoother.periods;
opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods;
+opts_simul.max_check_ahead_periods = options_.occbin.smoother.max_check_ahead_periods;
opts_simul.periodic_solution = options_.occbin.smoother.periodic_solution;
opts_simul.full_output = options_.occbin.smoother.full_output;
opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only;
@@ -254,7 +255,7 @@ while is_changed && maxiter>iter && ~is_periodic
err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
end
- if occbin_smoother_debug
+ if occbin_smoother_debug || is_periodic
regime_ = cell(0);
regime_new = regime_;
start_ = regime_;
@@ -332,18 +333,26 @@ regime_history0(max(iter+1,1),:) = regime_history;
oo_.occbin.smoother.regime_history=regime_history0(end,:);
oo_.occbin.smoother.regime_history_iter=regime_history0;
if occbin_smoother_debug
-save('info1','regime_history0')
+ save('info1','regime_history0')
end
if (maxiter==iter && is_changed) || is_periodic
- disp(['Occbin smoother did not converge.'])
+ disp('occbin.DSGE_smoother: smoother did not converge.')
+ fprintf('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.\n')
if is_periodic
- disp(['Occbin smoother algo loops between two solutions.'])
+ oo_.occbin.smoother.error_flag=0;
+ fprintf('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.\n')
+ fprintf('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.\n')
+ else
+ fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n')
+ oo_.occbin.smoother=[];
+ oo_.occbin.smoother.error_flag=1;
end
else
- disp(['Occbin smoother converged.'])
+ disp('occbin.DSGE_smoother: smoother converged.')
+ oo_.occbin.smoother.error_flag=0;
if occbin_smoother_fast && is_changed_start
- disp('WARNING: fast algo is used, regime(s) duration(s) was not forced to converge')
+ disp('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge')
end
end
if (~is_changed || occbin_smoother_debug) && nargin==12
diff --git a/matlab/+occbin/IVF_core.m b/matlab/+occbin/IVF_core.m
index 31e56e3c0f712d57b26e062f1a5412bf538b4aba..1ffe0a153515810cdc992906a6b0cae5fc17b1de 100644
--- a/matlab/+occbin/IVF_core.m
+++ b/matlab/+occbin/IVF_core.m
@@ -60,6 +60,7 @@ opts_simul.maxit = options_.occbin.likelihood.maxit;
opts_simul.waitbar = false;
opts_simul.periods = options_.occbin.likelihood.periods;
opts_simul.check_ahead_periods = options_.occbin.likelihood.check_ahead_periods;
+opts_simul.max_check_ahead_periods = options_.occbin.likelihood.max_check_ahead_periods;
opts_simul.periodic_solution = options_.occbin.likelihood.periodic_solution;
opts_simul.restrict_state_space = options_.occbin.likelihood.restrict_state_space;
opts_simul.piecewise_only = 1;
diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index ea02b1a1f2bc83c0d996f6056eff0fe2291c0c53..68d90928c792d8a40b942f3c855dde8a7ba5c615 100644
--- a/matlab/+occbin/set_default_options.m
+++ b/matlab/+occbin/set_default_options.m
@@ -76,6 +76,7 @@ if ismember(flag,{'likelihood','all'})
options_occbin_.likelihood.IVF_shock_observable_mapping = [];
options_occbin_.likelihood.maxit = 30; % this is for occbin solver algo
options_occbin_.likelihood.max_number_of_iterations = 10; % this is for occbin_kalman_update loop
+ options_occbin_.likelihood.max_check_ahead_periods=inf;
options_occbin_.likelihood.periods = 100;
options_occbin_.likelihood.check_ahead_periods=200;
options_occbin_.likelihood.periodic_solution=false;
@@ -170,11 +171,12 @@ if ismember(flag,{'simul','all'})
options_occbin_.simul.init_binding_indicator=false(0);
options_occbin_.simul.exo_pos=1:M_.exo_nbr;
options_occbin_.simul.maxit=30;
- options_occbin_.simul.max_periods=inf;
+ options_occbin_.simul.max_check_ahead_periods=inf;
options_occbin_.simul.periods=100;
options_occbin_.simul.check_ahead_periods=200;
options_occbin_.simul.periodic_solution=false;
options_occbin_.simul.piecewise_only = false;
+ options_occbin_.simul.reset_check_ahead_periods_in_new_period = false;
options_occbin_.simul.reset_regime_in_new_period = false;
options_occbin_.simul.restrict_state_space=false;
options_occbin_.simul.SHOCKS=zeros(1,M_.exo_nbr);
@@ -193,6 +195,7 @@ if ismember(flag,{'smoother','all'})
options_occbin_.smoother.inversion_filter = false;
options_occbin_.smoother.linear_smoother = true;
options_occbin_.smoother.maxit = 30; % this is for occbin solver algo
+ options_occbin_.smoother.max_check_ahead_periods=inf;
options_occbin_.smoother.max_number_of_iterations = 10; % this is for smoother loop
options_occbin_.smoother.periods = 100;
options_occbin_.smoother.check_ahead_periods=200;
diff --git a/matlab/+occbin/setup.m b/matlab/+occbin/setup.m
index 77c984c22f7435ef1db9b9145987aee5d175dd54..67db5bdb14dd6bbc95860f53d769df5820117448 100644
--- a/matlab/+occbin/setup.m
+++ b/matlab/+occbin/setup.m
@@ -2,12 +2,12 @@ function [M_, options_] = setup(M_,options_, options_occbin_)
% function [M_, options_] = setup(M_, options_, options_occbin_)
% Sets up run of Occbin: creates shock matrix, sets options
%
-% INPUT:
+% INPUT:
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure containing the options
% - options_occbin_ [structure] Matlab's structure containing Occbin options
%
-% OUTPUT:
+% OUTPUT:
% - M_ [structure] Matlab's structure describing the model
% - options_occbin_ [structure] Matlab's structure containing Occbin options
@@ -55,7 +55,7 @@ if isfield(M_,'surprise_shocks') && ~isempty(M_.surprise_shocks)
for ii = 1:length(M_.surprise_shocks)
ivar = M_.surprise_shocks(ii).exo_id;
temp(M_.surprise_shocks(ii).periods,ivar) = M_.surprise_shocks(ii).value;
- end
+ end
shock_index=~all(temp==0,1);
options_.occbin.simul.SHOCKS=temp(:,shock_index);
options_.occbin.simul.exo_pos=find(shock_index);
diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index 0801f8c1ebfab61a3b8e4658b30fee5bd22454a2..f5f0c1446c79ff6edf9b1c829e8b98c90f2216ca 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -125,21 +125,29 @@ 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
+ if binding_indicator(end) && nperiods_0<opts_simul_.max_check_ahead_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_check_ahead_periods
+ % enforce endogenously increased nperiods to not violate max_check_ahead_periods
+ binding_indicator = binding_indicator(1:opts_simul_.max_check_ahead_periods+1);
+ binding_indicator(end)=false;
end
if length(binding_indicator)<(nperiods_0 + 1)
binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)];
@@ -164,7 +172,13 @@ for shock_period = 1:n_shocks_periods
% get the hypothesized piece wise linear solution
if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:))
if iter==1 && opts_simul_.reset_regime_in_new_period
- binding_indicator=false(size(binding_indicator));
+ if opts_simul_.reset_check_ahead_periods_in_new_period
+ % I re-set check ahead periods to initial value, when in previous period it was endogenously increased
+ nperiods_0 = max(opts_simul_.check_ahead_periods,n_periods-n_shocks_periods);
+ binding_indicator = false(nperiods_0+1,1);
+ else
+ binding_indicator=false(size(binding_indicator));
+ end
binding_indicator_history{iter}=binding_indicator;
% analyse violvec and isolate contiguous periods in the other regime.
[regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
@@ -178,42 +192,81 @@ for shock_period = 1:n_shocks_periods
[binding, relax, err]=feval([M_.fname,'.occbin_difference'],zdatalinear_+repmat(dr_base.ys',size(zdatalinear_,1),1),M_.params,dr_base.ys);
+ if ~isinf(opts_simul_.max_check_ahead_periods) && opts_simul_.max_check_ahead_periods<length(binding_indicator)
+ end_periods = opts_simul_.max_check_ahead_periods;
+ last_indicator = false(length(binding_indicator)-end_periods,1);
+ else
+ end_periods = length(binding_indicator);
+ last_indicator = false(0);
+ end
% check if changes to the hypothesis of the duration for each
% regime
- if any(binding.constraint_1 & ~binding_indicator) || any(relax.constraint_1 & binding_indicator)
- err_viol = err.binding_constraint_1(binding.constraint_1 & ~binding_indicator);
- err_relax = err.relax_constraint_1(relax.constraint_1 & binding_indicator);
+ if any(binding.constraint_1(1:end_periods) & ~binding_indicator(1:end_periods)) || any(relax.constraint_1(1:end_periods) & binding_indicator(1:end_periods))
+
+ err_viol = err.binding_constraint_1(binding.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]));
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;
end
+ % if max_check_ahead_periods<inf, enforce unconstrained regime for period larger than max_check_ahead_periods
+ binding_constraint_new=[binding.constraint_1(1:end_periods);last_indicator];
+ relaxed_constraint_new = [relax.constraint_1(1:end_periods);not(last_indicator)];
if curb_retrench % apply Gauss-Seidel idea of slowing down the change in the guess
% for the constraint -- only relax one
% period at a time starting from the last
% one when each of the constraints is true.
retrench = false(numel(binding_indicator),1);
- max_relax_constraint_1=find(relax.constraint_1 & binding_indicator,1,'last');
- if ~isempty(max_relax_constraint_1) && find(relax.constraint_1,1,'last')>=find(binding_indicator,1,'last')
+ max_relax_constraint_1=find(relax.constraint_1(1:end_periods) & binding_indicator(1:end_periods),1,'last');
+ if ~isempty(max_relax_constraint_1) && find(relax.constraint_1(1:end_periods),1,'last')>=find(binding_indicator(1:end_periods),1,'last')
retrench(max_relax_constraint_1) = true;
end
- binding_indicator = (binding_indicator | binding.constraint_1) & ~ retrench;
+ binding_indicator = (binding_indicator | binding_constraint_new) & ~ retrench;
else
- binding_indicator= (binding_indicator | binding.constraint_1) & ~(binding_indicator & relax.constraint_1);
+ binding_indicator = (binding_indicator | binding_constraint_new) & ~(binding_indicator & relaxed_constraint_new);
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
@@ -265,9 +318,14 @@ for shock_period = 1:n_shocks_periods
return
end
else
+ if is_periodic_loop
+ disp_verbose('Did not converge -- infinite loop of guess 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
- error_flag = 311;
return
end
else
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index 234aaa78b510dfc93cc14d6b6c4b1532fc5c621c..80744ac50943fa24b2e3ee285cb7ad3d3c587152 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -133,23 +133,35 @@ 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
+ if any(binding_indicator(end,:)) && nperiods_0<opts_simul_.max_check_ahead_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_check_ahead_periods
+ % enforce endogenously increased nperiods to not violate max_check_ahead_periods
+ binding_indicator = binding_indicator(1:opts_simul_.max_check_ahead_periods+1,:);
+ binding_indicator(end,:)=false(1,2);
end
if size(binding_indicator,1)<(nperiods_0 + 1)
+ % to ensure the simulation is run for the required nperiods
+ % even beyond max_check_ahead_periods: the latter controls check ahead periods
+ % and NOT how many periods we simulate after we are back to
+ % unconstrained regime (nperiods_0)
binding_indicator=[binding_indicator; false(nperiods_0 + 1-size(binding_indicator,1),2)];
end
binding_indicator_history{iter}=binding_indicator;
@@ -178,7 +190,13 @@ for shock_period = 1:n_shocks_periods
regime_history(shock_period).regimestart2 = regime_start_2;
if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:)) % first period or shock happening
if iter==1 && opts_simul_.reset_regime_in_new_period
- binding_indicator=false(size(binding_indicator));
+ if opts_simul_.reset_check_ahead_periods_in_new_period
+ % I re-set check ahead periods to initial value, when in previous period it was endogenously increased
+ nperiods_0 = max(opts_simul_.check_ahead_periods,n_periods-n_shocks_periods);
+ binding_indicator = false(nperiods_0+1,2);
+ else
+ binding_indicator=false(size(binding_indicator));
+ end
binding_indicator_history{iter}=binding_indicator;
% analyse violvec and isolate contiguous periods in the other regime.
[regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
@@ -195,53 +213,91 @@ for shock_period = 1:n_shocks_periods
data.exo_pos,data.shocks_sequence(shock_period,:),endo_init, update_flag);
[binding, relax, err]=feval([M_.fname,'.occbin_difference'],zdatalinear_+repmat(dr.ys',size(zdatalinear_,1),1),M_.params,dr.ys);
- binding_constraint_new=[binding.constraint_1;binding.constraint_2];
- relaxed_constraint_new = [relax.constraint_1;relax.constraint_2];
+
+ if ~isinf(opts_simul_.max_check_ahead_periods) && opts_simul_.max_check_ahead_periods<length(binding_indicator)
+ end_periods = opts_simul_.max_check_ahead_periods;
+ last_indicator = false(length(binding_indicator)-end_periods,1);
+ else
+ end_periods = length(binding_indicator);
+ last_indicator = false(0);
+ end
+ binding_constraint_new=[binding.constraint_1(1:end_periods); binding.constraint_2(1:end_periods)];
+ relaxed_constraint_new = [relax.constraint_1(1:end_periods); relax.constraint_2(1:end_periods)];
+ my_binding_indicator = binding_indicator(1:end_periods,:);
- err_binding_constraint_new = [err.binding_constraint_1; err.binding_constraint_2];
- err_relaxed_constraint_new = [err.relax_constraint_1; err.relax_constraint_2];
+ err_binding_constraint_new = [err.binding_constraint_1(1:end_periods); err.binding_constraint_2(1:end_periods)];
+ err_relaxed_constraint_new = [err.relax_constraint_1(1:end_periods); err.relax_constraint_2(1:end_periods)];
% check if changes_
- if any(binding_constraint_new & ~binding_indicator(:)) || any(relaxed_constraint_new & binding_indicator(:))
- err_violation = err_binding_constraint_new(binding_constraint_new & ~binding_indicator(:));
- err_relax = err_relaxed_constraint_new(relaxed_constraint_new & binding_indicator(:));
+ if any(binding_constraint_new & ~my_binding_indicator(:)) || any(relaxed_constraint_new & my_binding_indicator(:))
+ err_violation = err_binding_constraint_new(binding_constraint_new & ~my_binding_indicator(:));
+ err_relax = err_relaxed_constraint_new(relaxed_constraint_new & my_binding_indicator(:));
max_err(iter) = max(abs([err_violation;err_relax]));
regime_change_this_iteration = true;
+ regime_violates_constraint_in_expectation(iter) = any(binding_constraint_new & ~binding_indicator(:));
else
regime_change_this_iteration = false;
max_err(iter) = 0;
end
+ binding_constraint_new=[binding.constraint_1(1:end_periods);last_indicator; binding.constraint_2(1:end_periods);last_indicator];
+ relaxed_constraint_new = [relax.constraint_1(1:end_periods);not(last_indicator); relax.constraint_2(1:end_periods);not(last_indicator)];
if curb_retrench % apply Gauss-Seidel idea of slowing down the change in the guess
% for the constraint -- only relax one
% period at a time starting from the last
% one when each of the constraints is true.
retrench = false(numel(binding_indicator),1);
- max_relax_constraint_1=find(relax.constraint_1 & binding_indicator(:,1),1,'last');
- if ~isempty(max_relax_constraint_1) && find(relax.constraint_1,1,'last')>=find(binding_indicator(:,1),1,'last')
+ max_relax_constraint_1=find(relax.constraint_1(1:end_periods) & binding_indicator(1:end_periods,1),1,'last');
+ if ~isempty(max_relax_constraint_1) && find(relax.constraint_1(1:end_periods),1,'last')>=find(binding_indicator(1:end_periods,1),1,'last')
retrench(max_relax_constraint_1) = true;
end
- max_relax_constraint_2=find(relax.constraint_2 & binding_indicator(:,2),1,'last');
- if ~isempty(max_relax_constraint_2) && find(relax.constraint_2,1,'last')>=find(binding_indicator(:,2),1,'last')
+ max_relax_constraint_2=find(relax.constraint_2(1:end_periods) & binding_indicator(1:end_periods,2),1,'last');
+ if ~isempty(max_relax_constraint_2) && find(relax.constraint_2(1:end_periods),1,'last')>=find(binding_indicator(1:end_periods,2),1,'last')
retrench(max_relax_constraint_2+nperiods_0+1) = true;
end
binding_indicator = (binding_indicator(:) | binding_constraint_new) & ~ retrench;
else
- binding_indicator= (binding_indicator(:) | binding_constraint_new) & ~(binding_indicator(:) & relaxed_constraint_new);
+ binding_indicator = (binding_indicator(:) | binding_constraint_new) & ~(binding_indicator(:) & relaxed_constraint_new);
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
@@ -295,8 +351,13 @@ for shock_period = 1:n_shocks_periods
return;
end
else
- disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
- error_flag = 311;
+ if is_periodic_loop
+ disp_verbose('Did not converge -- infinite loop of guess 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/+occbin/solver.m b/matlab/+occbin/solver.m
index 6b08d8a1f15a9941e5bb6d0ab8253a721764ffc5..11ac70eafa860dc806753c2fa63a3e9bb7b2d1fc 100644
--- a/matlab/+occbin/solver.m
+++ b/matlab/+occbin/solver.m
@@ -68,6 +68,12 @@ else
oo_.dr=sto_dr;
end
+if options_.occbin.simul.check_ahead_periods>options_.occbin.simul.max_check_ahead_periods
+ options_.occbin.simul.check_ahead_periods=options_.occbin.simul.max_check_ahead_periods;
+ disp(['occbin::options::' simul '_check_ahead_periods cannot exceed ' simul '_max_check_ahead_periods'])
+ disp(['occbin::options::' simul '_check_ahead_periods is re-set to be equal to ' simul '_max_check_ahead_periods'])
+end
+
if M_.occbin.constraint_nbr==1
[out, ss, error_flag ] = occbin.solve_one_constraint(M_,oo_.dr,options_.occbin.simul,solve_dr);
elseif M_.occbin.constraint_nbr==2
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index c3e86a08e8a736af4ab298c534db6d44f7d3ecad..c31e1e6dd97d608e5b9d8add63e48b8310583f91 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -974,6 +974,7 @@ occbin_options.opts_simul.curb_retrench = DynareOptions.occbin.likelihood.curb_r
occbin_options.opts_simul.maxit = DynareOptions.occbin.likelihood.maxit;
occbin_options.opts_simul.periods = DynareOptions.occbin.likelihood.periods;
occbin_options.opts_simul.check_ahead_periods = DynareOptions.occbin.likelihood.check_ahead_periods;
+occbin_options.opts_simul.max_check_ahead_periods = DynareOptions.occbin.likelihood.max_check_ahead_periods;
occbin_options.opts_simul.periodic_solution = DynareOptions.occbin.likelihood.periodic_solution;
occbin_options.opts_simul.restrict_state_space = DynareOptions.occbin.likelihood.restrict_state_space;
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 5f2801a4ac91342027c0a1285ac8ea91577c1ad5..44a77ede2697d5b75f57b9cab722c6f4e678b942 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -185,11 +185,14 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && ~options_.
else
if options_.occbin.smoother.status
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
+ if oo_.occbin.smoother.error_flag==0
+ [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+ end
else
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
- end
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
end
+ end
if options_.forecast > 0
oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
end
@@ -597,8 +600,12 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
else
if options_.occbin.smoother.status
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = occbin.DSGE_smoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_,dataset_,dataset_info);
+ if oo_.occbin.smoother.error_flag==0
+ [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+ end
else
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
+ [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
end
[oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
end
diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m
index f1a1fe39473ecd505c354d2b904fbd40962d7207..fd94f8027bb5814b51397bd4e6c7ccf1f6d3460d 100644
--- a/matlab/evaluate_smoother.m
+++ b/matlab/evaluate_smoother.m
@@ -120,7 +120,9 @@ else
DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_);
end
if ~(options_.occbin.smoother.status && options_.occbin.smoother.inversion_filter)
- [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+ if ~options_.occbin.smoother.status || (options_.occbin.smoother.status && oo_.occbin.smoother.error_flag==0)
+ [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
+ end
else
[oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff);
end
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 1013184d69aeabf1c1434238082a4a718742bfee..72b95ac76f6f691d50c3e341797945960afbf04a 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 guess regimes';
case 320
message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.';
otherwise
diff --git a/matlab/plot_shock_decomposition.m b/matlab/plot_shock_decomposition.m
index ac8a8cbb1deb537175ad4417737f03f008e598a3..2bdb0c6ebc97fa7f945c52c4cb13e68766b2e437 100644
--- a/matlab/plot_shock_decomposition.m
+++ b/matlab/plot_shock_decomposition.m
@@ -470,13 +470,14 @@ switch type
q2a.aux.yss=steady_state_aux;
end
i_var0 = i_var;
+ steady_state_0 = steady_state;
[za, endo_names, endo_names_tex, steady_state, i_var, oo_] = ...
annualized_shock_decomposition(z,M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a);
if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp.use_shock_groups)
mygroup = options_.plot_shock_decomp.use_shock_groups;
options_.plot_shock_decomp.use_shock_groups='';
zafull = ...
- annualized_shock_decomposition(zfull(i_var0,:,:),M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a);
+ annualized_shock_decomposition(zfull(i_var0,:,:),M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state_0(i_var0),q2a);
options_.plot_shock_decomp.use_shock_groups = mygroup;
end
end
diff --git a/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc b/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
index 768dac606c5471ac7826bd1cb92d069703e21c60..9f504605efd4f75964e107a9d7d05db38afcdcdf 100644
--- a/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
+++ b/tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
@@ -30,7 +30,7 @@ steady;
% Option=2: random simulation
@#for option_val in [1, 2]
-
+
option=@{option_val};
%%%%%%%%%%%%%%%% Inputs stop here %%%%%%%%%%%%%%%%%%%%%
@@ -41,28 +41,28 @@ if option==1
periods 1:9, 10, 50, 90, 130, 131:169;
values -0.0001, -0.01,-0.02, 0.01, 0.02, 0;
end;
-
-
+
+
elseif option==2
- nperiods = 100;
+ nperiods = 40;
randn('seed',1);
- shockssequence = 1*randn(nperiods,1)*0.02 ;
-
+ shockssequence = 1*randn(nperiods,1)*0.02 ;
+
shocks(surprise,overwrite);
var erra;
- periods 1:100;
+ periods 1:40;
values (shockssequence);
end;
end
-% set inputs
+% set inputs
occbin_setup(smoother_debug);
occbin_solver(simul_debug,simul_periodic_solution,simul_periods=200,simul_maxit=300,simul_curb_retrench,simul_check_ahead_periods=300);
-occbin_write_regimes(filename='test',periods=[1:100]);
+occbin_write_regimes(filename='test',periods=[1:40]);
-%% Modify to plot IRFs
+%% Modify to plot IRFs
titlelist = char('c','lambdak','k','i','a');
percent = 'Percent';
value = 'value';
@@ -85,7 +85,7 @@ occbin_graph(noconstant) c erra lambdak k i a k;
end;
save('datasim.mat','c');
varobs c;
-
+
occbin_solver(simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200);
occbin_setup(smoother_periods=200,smoother_maxit=200,smoother_curb_retrench,smoother_check_ahead_periods=200);
calib_smoother(datafile=datasim);
@@ -110,4 +110,4 @@ occbin_graph(noconstant) c erra lambdak k i a k;
line2=100*[oo_.occbin.linear_smoother.SmoothedVariables.c-oo_.occbin.endo_ss.c,oo_.occbin.linear_smoother.SmoothedVariables.lambdak/100,oo_.occbin.linear_smoother.SmoothedVariables.k-oo_.occbin.endo_ss.k,oo_.occbin.linear_smoother.SmoothedVariables.i-oo_.occbin.endo_ss.i,oo_.occbin.linear_smoother.SmoothedVariables.a-oo_.occbin.endo_ss.a, oo_.occbin.linear_smoother.SmoothedShocks.erra/100];
occbin.make_chart(titlelist,legendlist,figtitle,ylabels,cat(3,line1,line2));
@#endif
-@#endfor
\ No newline at end of file
+@#endfor