Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
Loading items
Show changes
Showing
with 632 additions and 187 deletions
...@@ -105,7 +105,7 @@ end ...@@ -105,7 +105,7 @@ end
SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods); SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods);
SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods); SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods);
SS_out.C = nan(DM.n_vars,n_shocks_periods); SS_out.C = nan(DM.n_vars,n_shocks_periods);
if ~exist('regime_history_','var') || isempty(regime_history_guess) if isempty(regime_history_guess)
regime_history = struct(); regime_history = struct();
guess_history = false; guess_history = false;
else else
...@@ -125,28 +125,34 @@ for shock_period = 1:n_shocks_periods ...@@ -125,28 +125,34 @@ 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_check_ahead_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_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 end
if length(binding_indicator)<(nperiods_0 + 1) if length(binding_indicator)<(nperiods_0 + 1)
binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)]; binding_indicator=[binding_indicator; false(nperiods_0 + 1-length(binding_indicator),1)];
end end
binding_indicator_history{iter}=binding_indicator;
if iter==1 && guess_history_it if iter==1 && guess_history_it
regime = regime_history_guess(shock_period).regime; regime = regime_history_guess(shock_period).regime;
regime_start = regime_history_guess(shock_period).regimestart; regime_start = regime_history_guess(shock_period).regimestart;
...@@ -156,6 +162,7 @@ for shock_period = 1:n_shocks_periods ...@@ -156,6 +162,7 @@ for shock_period = 1:n_shocks_periods
end end
nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required
end end
binding_indicator_history{iter}=binding_indicator;
% analyze when each regime starts based on current guess % analyze when each regime starts based on current guess
[regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug); [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
regime_history(shock_period).regime = regime; regime_history(shock_period).regime = regime;
...@@ -164,7 +171,13 @@ for shock_period = 1:n_shocks_periods ...@@ -164,7 +171,13 @@ for shock_period = 1:n_shocks_periods
% get the hypothesized piece wise linear solution % get the hypothesized piece wise linear solution
if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:)) if shock_period==1 || shock_period>1 && any(data.shocks_sequence(shock_period,:))
if iter==1 && opts_simul_.reset_regime_in_new_period if iter==1 && opts_simul_.reset_regime_in_new_period
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)); binding_indicator=false(size(binding_indicator));
end
binding_indicator_history{iter}=binding_indicator; binding_indicator_history{iter}=binding_indicator;
% analyse violvec and isolate contiguous periods in the other regime. % analyse violvec and isolate contiguous periods in the other regime.
[regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug); [regime, regime_start, error_code_period]=occbin.map_regime(binding_indicator,opts_simul_.debug);
...@@ -178,42 +191,81 @@ for shock_period = 1:n_shocks_periods ...@@ -178,42 +191,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); [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 % check if changes to the hypothesis of the duration for each
% regime % regime
if any(binding.constraint_1 & ~binding_indicator) || any(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 & ~binding_indicator);
err_relax = err.relax_constraint_1(relax.constraint_1 & binding_indicator); 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])); 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;
end 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 if curb_retrench % apply Gauss-Seidel idea of slowing down the change in the guess
% for the constraint -- only relax one % for the constraint -- only relax one
% period at a time starting from the last % period at a time starting from the last
% one when each of the constraints is true. % one when each of the constraints is true.
retrench = false(numel(binding_indicator),1); retrench = false(numel(binding_indicator),1);
max_relax_constraint_1=find(relax.constraint_1 & 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,'last')>=find(binding_indicator,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; retrench(max_relax_constraint_1) = true;
end end
binding_indicator = (binding_indicator | binding.constraint_1) & ~ retrench; binding_indicator = (binding_indicator | binding_constraint_new) & ~ retrench;
else 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 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
...@@ -264,10 +316,15 @@ for shock_period = 1:n_shocks_periods ...@@ -264,10 +316,15 @@ for shock_period = 1:n_shocks_periods
error_flag = 310; error_flag = 310;
return return
end end
else
if is_periodic_loop
disp_verbose('Did not converge -- infinite loop of guess regimes.',opts_simul_.debug)
error_flag = 313;
else 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
......
...@@ -115,7 +115,7 @@ end ...@@ -115,7 +115,7 @@ end
SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods); SS_out.T = NaN(DM.n_vars,DM.n_vars,n_shocks_periods);
SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods); SS_out.R = NaN(DM.n_vars,DM.n_exo,n_shocks_periods);
SS_out.C = NaN(DM.n_vars,n_shocks_periods); SS_out.C = NaN(DM.n_vars,n_shocks_periods);
if ~exist('regime_history_','var') || isempty(regime_history_guess) if isempty(regime_history_guess)
regime_history = struct(); regime_history = struct();
guess_history = false; guess_history = false;
else else
...@@ -133,26 +133,37 @@ for shock_period = 1:n_shocks_periods ...@@ -133,26 +133,37 @@ 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_check_ahead_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_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 end
if size(binding_indicator,1)<(nperiods_0 + 1) 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)]; binding_indicator=[binding_indicator; false(nperiods_0 + 1-size(binding_indicator,1),2)];
end end
binding_indicator_history{iter}=binding_indicator;
if iter==1 && guess_history_it if iter==1 && guess_history_it
regime_1 = regime_history_guess(shock_period).regime1; regime_1 = regime_history_guess(shock_period).regime1;
...@@ -169,6 +180,8 @@ for shock_period = 1:n_shocks_periods ...@@ -169,6 +180,8 @@ for shock_period = 1:n_shocks_periods
end end
nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required nperiods_0 = size(binding_indicator,1)-1; %if history is present, update may be required
end end
binding_indicator_history{iter}=binding_indicator;
% analyse violvec and isolate contiguous periods in the other regime. % 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); [regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
regime_history(shock_period).regime1 = regime_1; regime_history(shock_period).regime1 = regime_1;
...@@ -178,7 +191,13 @@ for shock_period = 1:n_shocks_periods ...@@ -178,7 +191,13 @@ for shock_period = 1:n_shocks_periods
regime_history(shock_period).regimestart2 = regime_start_2; 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 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 if iter==1 && opts_simul_.reset_regime_in_new_period
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)); binding_indicator=false(size(binding_indicator));
end
binding_indicator_history{iter}=binding_indicator; binding_indicator_history{iter}=binding_indicator;
% analyse violvec and isolate contiguous periods in the other regime. % 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); [regime_1, regime_start_1, error_code_period(1)]=occbin.map_regime(binding_indicator(:,1),opts_simul_.debug);
...@@ -195,34 +214,46 @@ for shock_period = 1:n_shocks_periods ...@@ -195,34 +214,46 @@ for shock_period = 1:n_shocks_periods
data.exo_pos,data.shocks_sequence(shock_period,:),endo_init, update_flag); 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, 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];
err_binding_constraint_new = [err.binding_constraint_1; err.binding_constraint_2]; if ~isinf(opts_simul_.max_check_ahead_periods) && opts_simul_.max_check_ahead_periods<length(binding_indicator)
err_relaxed_constraint_new = [err.relax_constraint_1; err.relax_constraint_2]; 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(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_ % check if changes_
if any(binding_constraint_new & ~binding_indicator(:)) || any(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 & ~binding_indicator(:)); err_violation = err_binding_constraint_new(binding_constraint_new & ~my_binding_indicator(:));
err_relax = err_relaxed_constraint_new(relaxed_constraint_new & binding_indicator(:)); err_relax = err_relaxed_constraint_new(relaxed_constraint_new & my_binding_indicator(:));
max_err(iter) = max(abs([err_violation;err_relax])); max_err(iter) = max(abs([err_violation;err_relax]));
regime_change_this_iteration = true; regime_change_this_iteration = true;
regime_violates_constraint_in_expectation(iter) = any(binding_constraint_new & ~binding_indicator(:));
else else
regime_change_this_iteration = false; regime_change_this_iteration = false;
max_err(iter) = 0; max_err(iter) = 0;
end 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 if curb_retrench % apply Gauss-Seidel idea of slowing down the change in the guess
% for the constraint -- only relax one % for the constraint -- only relax one
% period at a time starting from the last % period at a time starting from the last
% one when each of the constraints is true. % one when each of the constraints is true.
retrench = false(numel(binding_indicator),1); retrench = false(numel(binding_indicator),1);
max_relax_constraint_1=find(relax.constraint_1 & 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,'last')>=find(binding_indicator(:,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; retrench(max_relax_constraint_1) = true;
end end
max_relax_constraint_2=find(relax.constraint_2 & 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,'last')>=find(binding_indicator(:,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; retrench(max_relax_constraint_2+nperiods_0+1) = true;
end end
binding_indicator = (binding_indicator(:) | binding_constraint_new) & ~ retrench; binding_indicator = (binding_indicator(:) | binding_constraint_new) & ~ retrench;
...@@ -231,17 +262,43 @@ for shock_period = 1:n_shocks_periods ...@@ -231,17 +262,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
...@@ -294,9 +351,14 @@ for shock_period = 1:n_shocks_periods ...@@ -294,9 +351,14 @@ for shock_period = 1:n_shocks_periods
if opts_simul_.waitbar; dyn_waitbar_close(hh); end if opts_simul_.waitbar; dyn_waitbar_close(hh); end
return; return;
end end
else
if is_periodic_loop
disp_verbose('Did not converge -- infinite loop of guess regimes.',opts_simul_.debug)
error_flag = 313;
else 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
......
function [oo_, out, ss] = solver(M_,oo_,options_) function [oo_, out, ss] = solver(M_,oo_,options_)
% function [oo_, out, ss] = solver(M_,oo_,options_,opts_simul) % function [oo_, out, ss] = solver(M_,oo_,options_)
% Solves the model with an OBC and produces simulations/IRFs % Solves the model with an OBC and produces simulations/IRFs
% %
% INPUT: % INPUT:
% - opts_simul [structure] Occbin simulation options
% - M_ [structure] Matlab's structure describing the model % - M_ [structure] Matlab's structure describing the model
% - oo_ [structure] Matlab's structure containing the results
% - options_ [structure] Matlab's structure containing the options % - options_ [structure] Matlab's structure containing the options
% %
% OUTPUT: % OUTPUT:
...@@ -51,7 +51,19 @@ else ...@@ -51,7 +51,19 @@ else
end end
if solve_dr if solve_dr
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
if options_.order>1
options_.order = 1;
end
[dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_); [dr,error_flag,M_,oo_] = resol(0,M_,options_,oo_);
out.error_flag=error_flag;
if error_flag
print_info(error_flag, options_.noprint, options_)
return;
end
oo_.dr = dr; oo_.dr = dr;
sto_dr=dr; sto_dr=dr;
sto_M=M_; sto_M=M_;
...@@ -59,6 +71,12 @@ else ...@@ -59,6 +71,12 @@ else
oo_.dr=sto_dr; oo_.dr=sto_dr;
end 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 if M_.occbin.constraint_nbr==1
[out, ss, error_flag ] = occbin.solve_one_constraint(M_,oo_.dr,options_.occbin.simul,solve_dr); [out, ss, error_flag ] = occbin.solve_one_constraint(M_,oo_.dr,options_.occbin.simul,solve_dr);
elseif M_.occbin.constraint_nbr==2 elseif M_.occbin.constraint_nbr==2
...@@ -68,16 +86,22 @@ end ...@@ -68,16 +86,22 @@ end
out.error_flag=error_flag; out.error_flag=error_flag;
if error_flag if error_flag
print_info(error_flag, options_.noprint, options_) print_info(error_flag, options_.noprint, options_)
out=[];
return; return;
end end
% add back steady state % add back steady state
if ~options_.occbin.simul.piecewise_only if ~options_.occbin.simul.piecewise_only
if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
out.linear = bsxfun(@plus, out.linear, out.ys');
else
out.linear = out.linear + out.ys'; out.linear = out.linear + out.ys';
end end
end
if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
out.piecewise = bsxfun(@plus, out.piecewise, out.ys');
else
out.piecewise = out.piecewise + out.ys'; out.piecewise = out.piecewise + out.ys';
out.exo_simul = options_.occbin.simul.SHOCKS; end
out.exo_pos = options_.occbin.simul.exo_pos; out.exo_pos = options_.occbin.simul.exo_pos;
oo_.occbin=out; oo_.occbin.simul=out;
...@@ -29,10 +29,10 @@ function oo_=unpack_simulations(M_,oo_,options_) ...@@ -29,10 +29,10 @@ function oo_=unpack_simulations(M_,oo_,options_)
for i=1:M_.endo_nbr for i=1:M_.endo_nbr
% unpack the IRFs % unpack the IRFs
oo_.occbin.endo_linear.(M_.endo_names{i})= oo_.occbin.linear(:,i); oo_.occbin.endo_linear.(M_.endo_names{i})= oo_.occbin.simul.linear(:,i);
oo_.occbin.endo_piecewise.(M_.endo_names{i})=oo_.occbin.piecewise(:,i); oo_.occbin.endo_piecewise.(M_.endo_names{i})=oo_.occbin.simul.piecewise(:,i);
oo_.occbin.endo_ss.(M_.endo_names{i})=oo_.occbin.ys(i); oo_.occbin.endo_ss.(M_.endo_names{i})=oo_.occbin.simul.ys(i);
end end
for i=1:length(oo_.occbin.exo_pos) for i=1:length(oo_.occbin.simul.exo_pos)
oo_.occbin.exo.(M_.exo_names{i})=options_.occbin.simul.SHOCKS(:,i); oo_.occbin.exo.(M_.exo_names{i})=options_.occbin.simul.SHOCKS(:,i);
end end
\ No newline at end of file
function write_regimes_to_xls(regime_history,M_,options_) function write_regimes_to_xls(occbin_struct,M_,options_)
% function write_regimes_to_xls(regime_history,M_,options_) % function write_regimes_to_xls(occbin_struct,M_,options_)
% writes regime results to Excel-file % writes regime results to Excel-file
% %
% INPUTS % INPUTS
% - regime_history [struct] information on the regimes % - occbin_struct [struct] occbin structure containing information on the regimes
% - M_ [struct] Matlab's structure describing the model % - M_ [struct] Matlab's structure describing the model
% - options_ [struct] Matlab's structure describing the current options % - options_ [struct] Matlab's structure describing the current options
...@@ -26,6 +26,16 @@ function write_regimes_to_xls(regime_history,M_,options_) ...@@ -26,6 +26,16 @@ function write_regimes_to_xls(regime_history,M_,options_)
OutputDirectoryName = CheckPath('Output',M_.dname); OutputDirectoryName = CheckPath('Output',M_.dname);
if strcmpi(options_.occbin.write_regimes.type,'simul') || strcmpi(options_.occbin.write_regimes.type,'smoother')
if isfield(occbin_struct,options_.occbin.write_regimes.type) && isfield(occbin_struct.(options_.occbin.write_regimes.type),'regime_history')
regime_history=occbin_struct.(lower(options_.occbin.write_regimes.type)).regime_history;
else
error('write_regimes_to_xls: the required field does not exist');
end
else
error('write_regimes_to_xls: output type can only be simul or smoother.')
end
if isempty(options_.occbin.write_regimes.periods) if isempty(options_.occbin.write_regimes.periods)
T=1:length(regime_history); T=1:length(regime_history);
else else
...@@ -51,14 +61,30 @@ else ...@@ -51,14 +61,30 @@ else
xlsmat{tp,5}=int2str(regime_history(tp).regimestart2); xlsmat{tp,5}=int2str(regime_history(tp).regimestart2);
end end
end end
if ~ispc && ~isoctave && matlab_ver_less_than('9.0')
% On GNU/Linux and macOS, with MATLAB < R2016a, “writeable” can’t write Excel files
% (and “xlswrite” can’t either)
warning('This version of MATLAB is too old and cannot create Excel files. The Occbin regimes will rather be written to a CSV file.')
filename=[OutputDirectoryName filesep xls_filename '.csv'];
else
filename=[OutputDirectoryName filesep xls_filename '.xls']; filename=[OutputDirectoryName filesep xls_filename '.xls'];
if matlab_ver_less_than('9.3')
if exist(filename,'file')
delete(filename)
end end
else
if isfile(filename) if isfile(filename)
delete(filename) delete(filename)
end end
if isoctave
% “writetable” and “array2table” don’t exist under Octave
if isoctave && ~user_has_octave_forge_package('io')
error('The io package is required to write XLS files from Octave')
end end
xlswrite(filename, vertcat(Header, xlsmat));
elseif ~ispc && matlab_ver_less_than('9.0')
% Use a CSV file. See the comment above about filename.
% We don’t use xlswrite because its CSV fallback does not support cell-arrays.
writetable(array2table(xlsmat,'VariableNames',Header), filename);
else
writetable(array2table(xlsmat,'VariableNames',Header), filename, 'Sheet', 'Regimes'); writetable(array2table(xlsmat,'VariableNames',Header), filename, 'Sheet', 'Regimes');
end
...@@ -22,7 +22,7 @@ function iterative_ols(eqname, params, data, range) ...@@ -22,7 +22,7 @@ function iterative_ols(eqname, params, data, range)
% equation must have NaN values in the object. % equation must have NaN values in the object.
% [4] It is assumed that the residual is additive. % [4] It is assumed that the residual is additive.
% Copyright (C) 2018-2021 Dynare Team % Copyright (C) 2018-2022 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -121,42 +121,64 @@ end ...@@ -121,42 +121,64 @@ end
% Build PAC expectation matrix expression. % Build PAC expectation matrix expression.
dataForPACExpectation0 = dseries(); dataForPACExpectation0 = dseries();
listofvariables0 = {}; listofvariables0 = {};
isconstant = false;
if ~isempty(M_.pac.(pacmodl).equations.(eqtag).h0_param_indices) if ~isempty(M_.pac.(pacmodl).equations.(eqtag).h0_param_indices)
for i=1:length(M_.pac.(pacmodl).equations.(eqtag).h0_param_indices) for i=1:length(M_.pac.(pacmodl).equations.(eqtag).h0_param_indices)
match = regexp(rhs, sprintf('(?<var>((\\w*)|\\w*\\(-1\\)))\\*%s', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}), 'names'); match = regexp(rhs, sprintf('(?<var>((\\w*)|\\w*\\(-1\\)))\\*%s', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}), 'names');
if isempty(match) if isempty(match)
match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}), 'names'); match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}), 'names');
end end
if ~isempty(match)
if isempty(strfind(match.var, '(-1)')) if isempty(strfind(match.var, '(-1)'))
listofvariables0{i} = match.var; listofvariables0{end+1} = match.var;
dataForPACExpectation0 = [dataForPACExpectation0, data{listofvariables0{i}}]; dataForPACExpectation0 = [dataForPACExpectation0, data{listofvariables0{i}}];
else else
listofvariables0{i} = match.var(1:end-4); listofvariables0{end+1} = match.var(1:end-4);
dataForPACExpectation0 = [dataForPACExpectation0, data{match.var(1:end-4)}.lag(1)]; dataForPACExpectation0 = [dataForPACExpectation0, data{match.var(1:end-4)}.lag(1)];
end end
else
if strcmp(M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h0_param_indices(i)}, sprintf('h0_%s_eq0_constant', pacmodl))
isconstant = true;
end
end
end end
dataPAC0 = dataForPACExpectation0{listofvariables0{:}}(range).data; dataPAC0 = dataForPACExpectation0{listofvariables0{:}}(range).data;
if isconstant
dataPAC0 = [ones(rows(dataPAC0),1), dataPAC0];
end
else else
dataPAC0 = []; dataPAC0 = [];
end end
dataForPACExpectation1 = dseries(); dataForPACExpectation1 = dseries();
listofvariables1 = {}; listofvariables1 = {};
isconstant = false;
if ~isempty(M_.pac.(pacmodl).equations.(eqtag).h1_param_indices) if ~isempty(M_.pac.(pacmodl).equations.(eqtag).h1_param_indices)
for i=1:length(M_.pac.(pacmodl).equations.(eqtag).h1_param_indices) for i=1:length(M_.pac.(pacmodl).equations.(eqtag).h1_param_indices)
match = regexp(rhs, sprintf('(?<var>((\\w*)|(\\w*\\(-1\\))))\\*%s', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}), 'names'); match = regexp(rhs, sprintf('(?<var>((\\w*)|(\\w*\\(-1\\))))\\*%s', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}), 'names');
if isempty(match) if isempty(match)
match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}), 'names'); match = regexp(rhs, sprintf('%s\\*(?<var>((\\w*\\(-1\\))|(\\w*)))', M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}), 'names');
end end
match
if ~isempty(match)
if isempty(strfind(match.var, '(-1)')) if isempty(strfind(match.var, '(-1)'))
listofvariables1{i} = match.var; listofvariables1{end+1} = match.var;
dataForPACExpectation1 = [dataForPACExpectation1, data{listofvariables1{i}}]; dataForPACExpectation1 = [dataForPACExpectation1, data{listofvariables1{i}}];
else else
listofvariables1{i} = match.var(1:end-4); listofvariables1{end+1} = match.var(1:end-4);
dataForPACExpectation1 = [dataForPACExpectation1, data{match.var(1:end-4)}.lag(1)]; dataForPACExpectation1 = [dataForPACExpectation1, data{match.var(1:end-4)}.lag(1)];
end end
else
if strcmp(M_.param_names{M_.pac.(pacmodl).equations.(eqtag).h1_param_indices(i)}, sprintf('h1_%s_eq0_constant', pacmodl))
isconstant = true;
end end
end
end
listofvariables1
dataPAC1 = dataForPACExpectation1{listofvariables1{:}}(range).data; dataPAC1 = dataForPACExpectation1{listofvariables1{:}}(range).data;
if isconstant
dataPAC1 = [ones(rows(dataPAC1),1), dataPAC1];
end
else else
dataPAC1 = []; dataPAC1 = [];
end end
...@@ -303,6 +325,10 @@ end ...@@ -303,6 +325,10 @@ end
% Get indices in params0 for other parameters (optimizing agents share plus parameters related to exogenous variables). % Get indices in params0 for other parameters (optimizing agents share plus parameters related to exogenous variables).
[~, params_id_2] = setdiff(1:length(ipnames_), params_id_1); [~, params_id_2] = setdiff(1:length(ipnames_), params_id_1);
if isoctave && octave_ver_less_than('6')
% Under Octave < 6, setdiff() behaves as with the 'legacy' option (i.e. pre-R2012b behaviour)
params_id_2 = params_id_2';
end
% Get indices in params0 for the parameters associated to the exogenous variables. % Get indices in params0 for the parameters associated to the exogenous variables.
params_id_3 = setdiff(params_id_2, params_id_0); params_id_3 = setdiff(params_id_2, params_id_0);
......
...@@ -164,6 +164,13 @@ fprintf(fid, 'r = %s-(%s);\n', lhs, rhs); ...@@ -164,6 +164,13 @@ fprintf(fid, 'r = %s-(%s);\n', lhs, rhs);
fprintf(fid, 's = r''*r;\n'); fprintf(fid, 's = r''*r;\n');
fclose(fid); fclose(fid);
% Workaround for Octave bug https://savannah.gnu.org/bugs/?46282
% Octave will randomly fail to read the ssr_* file generated in the +folder
if isoctave
rename(['+' M_.fname], ['+' M_.fname '-tmp']);
rename(['+' M_.fname '-tmp'], ['+' M_.fname]);
end
% Copy (sub)sample data in a matrix. % Copy (sub)sample data in a matrix.
DATA = data([range(1)-1, range]).data; DATA = data([range(1)-1, range]).data;
......
...@@ -11,7 +11,7 @@ function parameters(pacname) ...@@ -11,7 +11,7 @@ function parameters(pacname)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2019 Dynare Team % Copyright © 2019-2021 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -69,8 +69,77 @@ for e=1:number_of_pac_eq ...@@ -69,8 +69,77 @@ for e=1:number_of_pac_eq
% Get PAC equation % Get PAC equation
pac_equation = equations.(eqtag); pac_equation = equations.(eqtag);
% Get Error correction and autoregressive parameters in PAC equation % Get Error correction and autoregressive parameters in PAC equation
a = NaN(1+pac_equation.max_lag, 1); params = NaN(2+pac_equation.max_lag, 1);
a(1) = M_.params(pac_equation.ec.params); params(1) = M_.params(pac_equation.ec.params);
a(1+(1:pac_equation.max_lag)) = M_.params(pac_equation.ar.params); params(1+(1:pac_equation.max_lag)) = M_.params(pac_equation.ar.params);
M_.params(pac_equation.mce.alpha) = a2alpha(a); params(end) = M_.params(pacmodel.discount_index);
[G, alpha, beta] = buildGmatrixWithAlphaAndBeta(params);
M_.params(pac_equation.mce.alpha) = flip(alpha);
if isfield(pacmodel, 'growth_neutrality_param_index')
if isfield(equations.(eqtag), 'non_optimizing_behaviour')
gamma = M_.params(equations.(eqtag).share_of_optimizing_agents_index);
else
gamma = 1.0;
end
A = [alpha; 1];
A_1 = polyval(A, 1.0);
A_b = polyval(A, beta);
m = length(alpha);
d = A_1*A_b*(iota(m, m)'*inv((eye(m)-G)*(eye(m)-G))*iota(m, m));
cc = 1/gamma-(sum(params(2:end-1))+d);
ll = 0;
if isfield(equations.(eqtag), 'optim_additive')
% Exogenous variables are present in the λ part (optimizing agents).
tmp0 = 0;
for i=1:length(equations.(eqtag).optim_additive.params)
if isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
tmp0 = tmp0 + equations.(eqtag).optim_additive.scaling_factor(i);
elseif ~isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
tmp0 = tmp0 + M_.params(equations.(eqtag).optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i);
elseif ~islogical(equations.(eqtag).optim_additive.bgp{i})
error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
end
end
cc = cc - tmp0;
end
if gamma<1
if isfield(equations.(eqtag), 'non_optimizing_behaviour') && isfield(equations.(eqtag).non_optimizing_behaviour, 'params')
% Exogenous variables are present in the 1-λ part (rule of thumb agents).
tmp0 = 0;
tmp1 = 0;
for i=1:length(equations.(eqtag).non_optimizing_behaviour.params)
if isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
tmp0 = tmp0 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
elseif ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
tmp0 = tmp0 + M_.params(equations.(eqtag).non_optimizing_behaviour.params(i))*equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.params(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
end
end
cc = cc - (1-gamma)*tmp0/gamma;
ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end
end
if isfield(equations.(eqtag), 'additive')
% Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
tmp0 = 0;
tmp1 = 0;
for i=1:length(equations.(eqtag).additive.params)
if isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
tmp0 = tmp0 + equations.(eqtag).additive.scaling_factor(i);
elseif ~isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
tmp0 = tmp0 + M_.params(equations.(eqtag).additive.params(i))*equations.(eqtag).additive.scaling_factor(i);
elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && isnan(equations.(eqtag).additive.params(i))
tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp{i};
elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && ~isnan(equations.(eqtag).additive.params(i))
tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.params(i)*equations.(eqtag).additive.bgp{i};
end
end
cc = cc - tmp0/gamma;
ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end
M_.params(pacmodel.growth_neutrality_param_index) = cc; % Multiplies the variable or expression provided though the growth option in command pac_model.
end
end end
\ No newline at end of file
...@@ -110,9 +110,15 @@ for e=1:number_of_pac_eq ...@@ -110,9 +110,15 @@ for e=1:number_of_pac_eq
if varmodel.nonstationary(id) if varmodel.nonstationary(id)
idns = id; idns = id;
ids = []; ids = [];
if varmodel.isconstant
idns = idns+1;
end
else else
idns = []; idns = [];
ids = id; ids = id;
if varmodel.isconstant
ids = ids+1;
end
end end
else else
% Trend component model is assumed. % Trend component model is assumed.
...@@ -127,6 +133,12 @@ for e=1:number_of_pac_eq ...@@ -127,6 +133,12 @@ for e=1:number_of_pac_eq
else else
growth_flag = false; growth_flag = false;
end end
% Do we have rule of thumb agents? γ is the share of optimizing agents.
if isfield(equations.(eqtag), 'non_optimizing_behaviour')
gamma = DynareModel.params(equations.(eqtag).share_of_optimizing_agents_index);
else
gamma = 1.0;
end
% Get h0 and h1 vectors (plus the parameter for the growth neutrality correction). % Get h0 and h1 vectors (plus the parameter for the growth neutrality correction).
if growth_flag if growth_flag
[h0, h1, growthneutrality] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, ids, idns, pacmodel.auxiliary_model_type); [h0, h1, growthneutrality] = hVectors([pacvalues; beta], varcalib.CompanionMatrix, ids, idns, pacmodel.auxiliary_model_type);
...@@ -169,56 +181,66 @@ for e=1:number_of_pac_eq ...@@ -169,56 +181,66 @@ for e=1:number_of_pac_eq
end end
end end
% Update the parameter related to the growth neutrality correction. % Update the parameter related to the growth neutrality correction.
if isfield(equations.(eqtag), 'non_optimizing_behaviour')
gamma = DynareModel.params(equations.(eqtag).share_of_optimizing_agents_index);
else
gamma = 1.0;
end
if growth_flag if growth_flag
% Growth neutrality as returned by hVector is valid iff % Growth neutrality as returned by hVector is valid iff
% there is no exogenous variables in the model and in the % there is no exogenous variables in the model and in the
% absence of non optimizing agents. % absence of non optimizing agents.
gg = -(growthneutrality-1); % Finite sum of autoregressive parameters + infinite sum of the coefficients in the PAC expectation term. gg = -(growthneutrality-1); % Finite sum of autoregressive parameters + infinite sum of the coefficients in the PAC expectation term.
cc = 1.0-gamma*gg; % First adjustment of the growth neutrality correction (should also be divided by gamma, done below at the end of this section). cc = 1.0-gg*gamma; % First adjustment of the growth neutrality correction (should also be divided by gamma, done below at the end of this section).
% We may have to further change the correction if we have nonzero mean exogenous variables. % We may have to further change the correction if we have nonzero mean exogenous variables.
ll = 0.0;
if isfield(equations.(eqtag), 'optim_additive') if isfield(equations.(eqtag), 'optim_additive')
% Exogenous variables are present in the λ part (optimizing agents). % Exogenous variables are present in the λ part (optimizing agents).
tmp = 0; tmp0 = 0;
for i=1:length(equations.(eqtag).optim_additive.params) for i=1:length(equations.(eqtag).optim_additive.params)
if isnan(equations.(eqtag).optim_additive.params(i)) && equations.(eqtag).optim_additive.bgp(i) if isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
tmp = tmp + equations.(eqtag).optim_additive.scaling_factor(i)*equations.(eqtag).optim_additive.bgp(i); tmp0 = tmp0 + equations.(eqtag).optim_additive.scaling_factor(i);
elseif ~isnan(equations.(eqtag).optim_additive.params(i)) && equations.(eqtag).optim_additive.bgp(i) elseif ~isnan(equations.(eqtag).optim_additive.params(i)) && islogical(equations.(eqtag).optim_additive.bgp{i}) && equations.(eqtag).optim_additive.bgp{i}
tmp = tmp + DynareModel.params(equations.(eqtag).optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i)*equations.(eqtag).optim_additive.bgp(i); tmp0 = tmp0 + DynareModel.params(equations.(eqtag).optim_additive.params(i))*equations.(eqtag).optim_additive.scaling_factor(i);
elseif ~islogical(equations.(eqtag).optim_additive.bgp{i})
error('It is not possible to provide a value for the mean of an exogenous variable appearing in the optimal part of the PAC equation.')
end end
end end
cc = cc - gamma*tmp; cc = cc - tmp0*gamma;
end end
if gamma<1 if gamma<1
if isfield(equations.(eqtag), 'non_optimizing_behaviour.params') if isfield(equations.(eqtag), 'non_optimizing_behaviour') && isfield(equations.(eqtag).non_optimizing_behaviour, 'params')
% Exogenous variables are present in the 1-λ part (rule of thumb agents). % Exogenous variables are present in the 1-λ part (rule of thumb agents).
tmp = 0; tmp0 = 0;
tmp1 = 0;
for i=1:length(equations.(eqtag).non_optimizing_behaviour.params) for i=1:length(equations.(eqtag).non_optimizing_behaviour.params)
if isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && equations.(eqtag).non_optimizing_behaviour.bgp(i) if isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
tmp = tmp + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp(i); tmp0 = tmp0 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
elseif ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && equations.(eqtag).non_optimizing_behaviour.bgp(i) elseif ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i)) && islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && equations.(eqtag).non_optimizing_behaviour.bgp{i}
tmp = tmp + DynareModel.params(equations.(eqtag).non_optimizing_behaviour.params(i))*equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp(i); tmp0 = tmp0 + DynareModel.params(equations.(eqtag).non_optimizing_behaviour.params(i))*equations.(eqtag).non_optimizing_behaviour.scaling_factor(i);
elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
elseif ~islogical(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && isnumeric(equations.(eqtag).non_optimizing_behaviour.bgp{i}) && ~isnan(equations.(eqtag).non_optimizing_behaviour.params(i))
tmp1 = tmp1 + equations.(eqtag).non_optimizing_behaviour.scaling_factor(i)*equations.(eqtag).non_optimizing_behaviour.params(i)*equations.(eqtag).non_optimizing_behaviour.bgp{i};
end end
end end
cc = cc - (1.0-gamma)*tmp; cc = cc - (1.0-gamma)*tmp0;
ll = -(1.0-gamma)*tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end end
end end
if isfield(equations.(eqtag), 'additive') if isfield(equations.(eqtag), 'additive')
% Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation. % Exogenous variables are present outside of the λ and (1-λ) parts (or we have exogenous variables in a "pure" PAC equation.
tmp = 0; tmp0 = 0;
tmp1 = 0;
for i=1:length(equations.(eqtag).additive.params) for i=1:length(equations.(eqtag).additive.params)
if isnan(equations.(eqtag).additive.params(i)) && equations.(eqtag).additive.bgp(i) if isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
tmp = tmp + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp(i); tmp0 = tmp0 + equations.(eqtag).additive.scaling_factor(i);
elseif ~isnan(equations.(eqtag).additive.params(i)) && equations.(eqtag).additive.bgp(i) elseif ~isnan(equations.(eqtag).additive.params(i)) && islogical(equations.(eqtag).additive.bgp{i}) && equations.(eqtag).additive.bgp{i}
tmp = tmp + DynareModel.params(equations.(eqtag).additive.params(i))*equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp(i); tmp0 = tmp0 + DynareModel.params(equations.(eqtag).additive.params(i))*equations.(eqtag).additive.scaling_factor(i);
elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && isnan(equations.(eqtag).additive.params(i))
tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.bgp{i};
elseif ~islogical(equations.(eqtag).additive.bgp{i}) && isnumeric(equations.(eqtag).additive.bgp{i}) && ~isnan(equations.(eqtag).additive.params(i))
tmp1 = tmp1 + equations.(eqtag).additive.scaling_factor(i)*equations.(eqtag).additive.params(i)*equations.(eqtag).additive.bgp{i};
end end
end end
cc = cc - tmp; cc = cc - tmp0;
ll = ll - tmp1/gamma; % TODO: ll should be added as a constant in the PAC equation (under the λ part) when unrolling pac_expectation.
end end
DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; DynareModel.params(pacmodel.growth_neutrality_param_index) = cc/gamma; % Multiplies the variable or expression provided though the growth option in command pac_model.
end end
end end
\ No newline at end of file
...@@ -8,7 +8,7 @@ function initialize(pacmodel) ...@@ -8,7 +8,7 @@ function initialize(pacmodel)
% OUTPUTS % OUTPUTS
% None % None
% Copyright (C) 2018-2019 Dynare Team % Copyright © 2018-2021 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -37,12 +37,12 @@ end ...@@ -37,12 +37,12 @@ end
% Append default bgp fields if needed. % Append default bgp fields if needed.
for i=1:rows(M_.pac.(pacmodel).tag_map) for i=1:rows(M_.pac.(pacmodel).tag_map)
if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'additive') if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'additive')
M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.bgp = false(1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.params)); M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).additive.params));
end end
if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'optim_additive') if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'optim_additive')
M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.bgp = false(1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.params)); M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).optim_additive.params));
end end
if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'non_optimizing_behaviour') if isfield(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}), 'non_optimizing_behaviour')
M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.bgp = false(1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.params)); M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.bgp = repmat({false}, 1, length(M_.pac.(pacmodel).equations.(M_.pac.(pacmodel).tag_map{i,2}).non_optimizing_behaviour.params));
end end
end end
\ No newline at end of file
...@@ -108,6 +108,10 @@ for i = 1:m ...@@ -108,6 +108,10 @@ for i = 1:m
end end
end end
if isfield(auxmodel, 'isconstant') && auxmodel.isconstant
variables_id_in_var = variables_id_in_var+1;
end
% Get the horizon parameter. % Get the horizon parameter.
horizon = varexpectationmodel.horizon; horizon = varexpectationmodel.horizon;
...@@ -175,8 +179,8 @@ else ...@@ -175,8 +179,8 @@ else
% Define the discounted companion matrix % Define the discounted companion matrix
DiscountedCompanionMatrix = discountfactor*CompanionMatrix; DiscountedCompanionMatrix = discountfactor*CompanionMatrix;
% First compute the parameters implied by the discounted sum from h=0 to h=horizon(1)-1 % First compute the parameters implied by the discounted sum from h=0 to h=horizon(1)-1
tmp1 = eye(n); tmp1 = zeros(n);
for h=1:horizon(1) for h=0:horizon(1)-1
tmp1 = tmp1 + mpower(DiscountedCompanionMatrix, h); tmp1 = tmp1 + mpower(DiscountedCompanionMatrix, h);
end end
tmp1 = alpha*tmp1; tmp1 = alpha*tmp1;
......
...@@ -189,6 +189,7 @@ elseif options_.lik_init == 2 % Old Diffuse Kalman filter ...@@ -189,6 +189,7 @@ elseif options_.lik_init == 2 % Old Diffuse Kalman filter
Pstar = options_.Harvey_scale_factor*eye(np); Pstar = options_.Harvey_scale_factor*eye(np);
Pinf = []; Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter elseif options_.lik_init == 3 % Diffuse Kalman filter
my_mf = mf;
if kalman_algo ~= 4 if kalman_algo ~= 4
kalman_algo = 3; kalman_algo = 3;
else else
...@@ -201,9 +202,10 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter ...@@ -201,9 +202,10 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
R = blkdiag(R,eye(vobs)); R = blkdiag(R,eye(vobs));
H = zeros(vobs,vobs); H = zeros(vobs,vobs);
Z = [Z, eye(vobs)]; Z = [Z, eye(vobs)];
my_mf = find(any(Z))';
end end
end end
[Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium); [Pstar,Pinf] = compute_Pinf_Pstar(my_mf,T,R,Q,options_.qz_criterium);
elseif options_.lik_init == 4 % Start from the solution of the Riccati equation. elseif options_.lik_init == 4 % Start from the solution of the Riccati equation.
Pstar = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H); Pstar = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
Pinf = []; Pinf = [];
...@@ -319,7 +321,7 @@ if kalman_algo == 2 || kalman_algo == 4 ...@@ -319,7 +321,7 @@ if kalman_algo == 2 || kalman_algo == 4
a_initial = zeros(np,1); a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_); a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=ST*a_initial; %set state prediction for first Kalman step; a_initial=ST*a_initial; %set state prediction for first Kalman step;
[alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ... [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
Z,R1,Q,diag(H), ... Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ... Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol, ... options_.nk,kalman_tol,diffuse_kalman_tol, ...
...@@ -363,6 +365,12 @@ else ...@@ -363,6 +365,12 @@ else
ic = oo_.dr.restrict_columns; ic = oo_.dr.restrict_columns;
end end
if isempty(options_.nk)
nk=1;
else
nk=options_.nk;
end
if options_.occbin.smoother.status if options_.occbin.smoother.status
% reconstruct occbin smoother % reconstruct occbin smoother
if length_varargin>0 if length_varargin>0
...@@ -386,9 +394,11 @@ else ...@@ -386,9 +394,11 @@ else
else else
opts_simul = options_.occbin.simul; opts_simul = options_.occbin.simul;
end end
% reconstruct smoothed variables
aaa=zeros(M_.endo_nbr,gend); aaa=zeros(M_.endo_nbr,gend);
aaa(oo_.dr.restrict_var_list,:)=alphahat; aaa(oo_.dr.restrict_var_list,:)=alphahat;
for k=2:gend iTx = zeros(size(TTx));
for k=1:gend
if isoccbin if isoccbin
A = TT(:,:,k); A = TT(:,:,k);
B = RR(:,:,k); B = RR(:,:,k);
...@@ -396,20 +406,53 @@ else ...@@ -396,20 +406,53 @@ else
else else
C=0; C=0;
end end
aaa(:,k) = C+A*aaa(:,k-1)+B*etahat(:,k); iT = pinv(TTx(:,:,k));
% store pinv
iTx(:,:,k) = iT;
Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
AS = Tstar*iT;
BS = Rstar-AS*RRx(:,:,k);
CS = Cstar-AS*CCx(:,k);
static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
static_var_list0 = static_var_list;
static_var_list0(static_var_list) = ilagged;
static_var_list(static_var_list) = ~ilagged;
aaa(static_var_list,k) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged);
if any(ilagged) && k>1
aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(ilagged);
end
end end
alphahat=aaa; alphahat=aaa;
aaa=zeros(M_.endo_nbr,gend);
% reconstruct updated variables
bbb=zeros(M_.endo_nbr,gend); bbb=zeros(M_.endo_nbr,gend);
bbb(oo_.dr.restrict_var_list,:)=ahat; bbb(oo_.dr.restrict_var_list,:)=ahat; % this is t|t
aaa(oo_.dr.restrict_var_list,:)=aahat; for k=1:gend
for k=d+2:gend
if isoccbin if isoccbin
A = TT(:,:,k); A = TT(:,:,k);
B = RR(:,:,k); B = RR(:,:,k);
C = CC(:,k); C = CC(:,k);
bbb(:,k) = C+A*aaa(:,k-1)+B*eehat(:,k); iT = iTx(:,:,k);
else Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
AS = Tstar*iT;
BS = Rstar-AS*RRx(:,:,k);
CS = Cstar-AS*CCx(:,k);
static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
static_var_list0 = static_var_list;
static_var_list0(static_var_list) = ilagged;
static_var_list(static_var_list) = ~ilagged;
bbb(static_var_list,k) = AS(~ilagged,:)*ahat(:,k)+BS(~ilagged,:)*eehat(:,k)+CS(~ilagged);
if any(ilagged) && k>d+1
bbb(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k)+Cstar(ilagged);
end
elseif k>d+1
opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench; opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench;
opts_simul.waitbar = options_.occbin.smoother.waitbar; opts_simul.waitbar = options_.occbin.smoother.waitbar;
opts_simul.maxit = options_.occbin.smoother.maxit; opts_simul.maxit = options_.occbin.smoother.maxit;
...@@ -417,23 +460,22 @@ else ...@@ -417,23 +460,22 @@ else
opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods; opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods;
opts_simul.full_output = options_.occbin.smoother.full_output; opts_simul.full_output = options_.occbin.smoother.full_output;
opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only; opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only;
opts_simul.SHOCKS = zeros(options_.nk,M_.exo_nbr); opts_simul.SHOCKS = zeros(nk,M_.exo_nbr);
opts_simul.SHOCKS(1,:) = eehat(:,k); opts_simul.SHOCKS(1,:) = eehat(:,k);
tmp=zeros(M_.endo_nbr,1); tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=aahat(:,k-1); tmp(oo_.dr.restrict_var_list,1)=aahat(:,k-1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1); opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
opts_simul.init_regime = []; %regimes_(k); opts_simul.init_regime = []; %regimes_(k);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out] = occbin.solver(M_,oo_,options_); [~, out] = occbin.solver(M_,oo_,options_);
% regime in out should be identical to regimes_(k-2) moved one % regime in out should be identical to regimes_(k-2) moved one
% period ahead (so if regimestart was [1 5] it should be [1 4] % period ahead (so if regimestart was [1 5] it should be [1 4]
% in out % in out
% end % end
bbb(oo_.dr.inv_order_var,k) = out.zpiece(1,:); bbb(oo_.dr.inv_order_var,k) = out.piecewise(1,:) - out.ys';
end end
end end
% do not overwrite accurate computations using reduced st. space
bbb(oo_.dr.restrict_var_list,:)=ahat;
ahat0=ahat; ahat0=ahat;
ahat=bbb; ahat=bbb;
if ~isempty(P) if ~isempty(P)
...@@ -444,13 +486,14 @@ else ...@@ -444,13 +486,14 @@ else
end end
if ~isempty(state_uncertainty) if ~isempty(state_uncertainty)
mm=size(T,1);
sstate_uncertainty=zeros(M_.endo_nbr,M_.endo_nbr,gend); sstate_uncertainty=zeros(M_.endo_nbr,M_.endo_nbr,gend);
sstate_uncertainty(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:)=state_uncertainty; sstate_uncertainty(oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:)=state_uncertainty(1:mm,1:mm,:);
state_uncertainty=sstate_uncertainty; state_uncertainty=sstate_uncertainty;
clear sstate_uncertainty clear sstate_uncertainty
end end
aaa = zeros(options_.nk,M_.endo_nbr,gend+options_.nk); aaa = zeros(nk,M_.endo_nbr,gend+nk);
aaa(:,oo_.dr.restrict_var_list,:)=aK; aaa(:,oo_.dr.restrict_var_list,:)=aK;
for k=2:gend+1 for k=2:gend+1
...@@ -461,25 +504,26 @@ else ...@@ -461,25 +504,26 @@ else
opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods; opts_simul.check_ahead_periods = options_.occbin.smoother.check_ahead_periods;
opts_simul.full_output = options_.occbin.smoother.full_output; opts_simul.full_output = options_.occbin.smoother.full_output;
opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only; opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only;
opts_simul.SHOCKS = zeros(options_.nk,M_.exo_nbr); opts_simul.SHOCKS = zeros(nk,M_.exo_nbr);
tmp=zeros(M_.endo_nbr,1); tmp=zeros(M_.endo_nbr,1);
tmp(oo_.dr.restrict_var_list,1)=ahat0(:,k-1); tmp(oo_.dr.restrict_var_list,1)=ahat0(:,k-1);
opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1); opts_simul.endo_init = tmp(oo_.dr.inv_order_var,1);
opts_simul.init_regime = []; %regimes_(k); opts_simul.init_regime = []; %regimes_(k);
opts_simul.waitbar=0;
options_.occbin.simul=opts_simul; options_.occbin.simul=opts_simul;
[~, out] = occbin.solver(M_,oo_,options_); [~, out] = occbin.solver(M_,oo_,options_);
% regime in out should be identical to regimes_(k-2) moved one % regime in out should be identical to regimes_(k-2) moved one
% period ahead (so if regimestart was [1 5] it should be [1 4] % period ahead (so if regimestart was [1 5] it should be [1 4]
% in out % in out
% end % end
for jnk=1:options_.nk for jnk=1:nk
aaa(jnk,oo_.dr.inv_order_var,k+jnk-1) = out.zpiece(jnk,:); aaa(jnk,oo_.dr.inv_order_var,k+jnk-1) = out.piecewise(jnk,:) - out.ys';
end end
end end
aK=aaa; aK=aaa;
if ~isempty(PK) if ~isempty(PK)
PP = zeros(options_.nk,M_.endo_nbr,M_.endo_nbr,gend+options_.nk); PP = zeros(nk,M_.endo_nbr,M_.endo_nbr,gend+nk);
PP(:,oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:) = PK; PP(:,oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:) = PK;
PK=PP; PK=PP;
clear PP clear PP
...@@ -525,17 +569,17 @@ else ...@@ -525,17 +569,17 @@ else
end end
ahat1=aaa; ahat1=aaa;
% reconstruct aK % reconstruct aK
aaa = zeros(options_.nk,M_.endo_nbr,gend+options_.nk); aaa = zeros(nk,M_.endo_nbr,gend+nk);
aaa(:,oo_.dr.restrict_var_list,:)=aK; aaa(:,oo_.dr.restrict_var_list,:)=aK;
for k=1:gend for k=1:gend
for jnk=1:options_.nk for jnk=1:nk
aaa(jnk,static_var_list,k+jnk) = C(~ilagged,:)*dynare_squeeze(aK(jnk,:,k+jnk)); aaa(jnk,static_var_list,k+jnk) = C(~ilagged,:)*dynare_squeeze(aK(jnk,:,k+jnk));
end end
end end
if any(ilagged) if any(ilagged)
for k=1:gend for k=1:gend
aaa(1,static_var_list0,k+1) = Tstar(ilagged,:)*ahat(:,k); aaa(1,static_var_list0,k+1) = Tstar(ilagged,:)*ahat(:,k);
for jnk=2:options_.nk for jnk=2:nk
aaa(jnk,static_var_list0,k+jnk) = Tstar(ilagged,:)*dynare_squeeze(aK(jnk-1,:,k+jnk-1)); aaa(jnk,static_var_list0,k+jnk) = Tstar(ilagged,:)*dynare_squeeze(aK(jnk-1,:,k+jnk-1));
end end
end end
...@@ -582,12 +626,12 @@ else ...@@ -582,12 +626,12 @@ else
% reconstruct PK % reconstruct PK
if ~isempty(PK) if ~isempty(PK)
PP = zeros(options_.nk,M_.endo_nbr,M_.endo_nbr,gend+options_.nk); PP = zeros(nk,M_.endo_nbr,M_.endo_nbr,gend+nk);
PP(:,oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:) = PK; PP(:,oo_.dr.restrict_var_list,oo_.dr.restrict_var_list,:) = PK;
if ~options_.heteroskedastic_filter if ~options_.heteroskedastic_filter
DQD=D(~ilagged,:)*Q*transpose(D(~ilagged,:))+C(~ilagged,:)*R*Q*transpose(D(~ilagged,:))+D(~ilagged,:)*Q*transpose(C(~ilagged,:)*R); DQD=D(~ilagged,:)*Q*transpose(D(~ilagged,:))+C(~ilagged,:)*R*Q*transpose(D(~ilagged,:))+D(~ilagged,:)*Q*transpose(C(~ilagged,:)*R);
DQR=D(~ilagged,:)*Q*transpose(R); DQR=D(~ilagged,:)*Q*transpose(R);
for f=1:options_.nk for f=1:nk
for k=1:gend for k=1:gend
PP(f,static_var_list,static_var_list,k+f)=C(~ilagged,:)*squeeze(PK(f,:,:,k+f))*C(~ilagged,:)'+DQD; PP(f,static_var_list,static_var_list,k+f)=C(~ilagged,:)*squeeze(PK(f,:,:,k+f))*C(~ilagged,:)'+DQD;
PP(f,static_var_list,oo_.dr.restrict_var_list,k+f)=C(~ilagged,:)*squeeze(PK(f,:,:,k+f))+DQR; PP(f,static_var_list,oo_.dr.restrict_var_list,k+f)=C(~ilagged,:)*squeeze(PK(f,:,:,k+f))+DQR;
......
...@@ -63,7 +63,7 @@ clear record; ...@@ -63,7 +63,7 @@ clear record;
header_width = row_header_width(M_, estim_params_, bayestopt_); header_width = row_header_width(M_, estim_params_, bayestopt_);
hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval']; hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval'];
tit2 = sprintf('%-*s %12s %12s %23s %8s %12s\n',header_width,' ','prior mean','post. mean',hpd_interval,'prior','pstdev'); tit2 = sprintf('%-*s %12s %12s %23s %8s %12s\n',header_width,' ','prior mean','post. mean',hpd_interval,'prior','pstdev');
pformat = '%-*s %12.3f % 12.4f %11.4f %11.4f %7s %12.4f'; pformat = '%-*s %12.3f % 12.4f %11.4f %11.4f %8s %12.4f';
skipline(2) skipline(2)
disp('ESTIMATION RESULTS') disp('ESTIMATION RESULTS')
......
...@@ -11,7 +11,7 @@ function WriteShockDecomp2Excel(z,shock_names,endo_names,i_var,initial_date,Dyna ...@@ -11,7 +11,7 @@ function WriteShockDecomp2Excel(z,shock_names,endo_names,i_var,initial_date,Dyna
% DynareModel [structure] Dynare model structure % DynareModel [structure] Dynare model structure
% DynareOptions [structure] Dynare options structure % DynareOptions [structure] Dynare options structure
% Copyright (C) 2016-2021 Dynare Team % Copyright (C) 2016-2022 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -124,7 +124,7 @@ for j=1:nvar ...@@ -124,7 +124,7 @@ for j=1:nvar
else else
writetable(cell2table(d0), [OutputDirectoryName,filesep,DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1 '.xls'], 'Sheet', endo_names{i_var(j)},'WriteVariableNames',false); writetable(cell2table(d0), [OutputDirectoryName,filesep,DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1 '.xls'], 'Sheet', endo_names{i_var(j)},'WriteVariableNames',false);
end end
warning on warning_config;
clear d0 clear d0
......
function mexpath = add_path_to_mex_files(dynareroot, modifypath) function mexpath = add_path_to_mex_files(dynareroot, modifypath)
% Copyright (C) 2015-2021 Dynare Team % Copyright © 2015-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -60,7 +60,7 @@ else ...@@ -60,7 +60,7 @@ else
end end
end end
else else
tmp = [dynareroot '../mex/matlab/win64-9.4-9.10/']; tmp = [dynareroot '../mex/matlab/win64-9.4-23.2/'];
if exist(tmp, 'dir') if exist(tmp, 'dir')
mexpath = tmp; mexpath = tmp;
if modifypath if modifypath
...@@ -80,7 +80,7 @@ else ...@@ -80,7 +80,7 @@ else
end end
end end
else else
tmp = [dynareroot '../mex/matlab/maci64-9.4-9.10']; tmp = [dynareroot '../mex/matlab/maci64-9.4-23.2/'];
if exist(tmp, 'dir') if exist(tmp, 'dir')
mexpath = tmp; mexpath = tmp;
if modifypath if modifypath
...@@ -89,6 +89,15 @@ else ...@@ -89,6 +89,15 @@ else
end end
end end
end end
if strcmp(computer, 'MACA64')
tmp = [dynareroot '../mex/matlab/maca64-23.2/'];
if exist(tmp, 'dir')
mexpath = tmp;
if modifypath
addpath(mexpath);
end
end
end
% Add generic MATLAB path (with higher priority than the previous ones) % Add generic MATLAB path (with higher priority than the previous ones)
if exist('mexpath') if exist('mexpath')
mexpath = { mexpath; [dynareroot '../mex/matlab/'] }; mexpath = { mexpath; [dynareroot '../mex/matlab/'] };
......
...@@ -21,7 +21,9 @@ function aggregate(ofile, dynopt, rootfolder, varargin) ...@@ -21,7 +21,9 @@ function aggregate(ofile, dynopt, rootfolder, varargin)
MAX_NUMBER_OF_ELEMENTS = 10000; MAX_NUMBER_OF_ELEMENTS = 10000;
if ~isoctave && matlab_ver_less_than('9.14') % Warning removed in R2023a
warning off MATLAB:subscripting:noSubscriptsSpecified warning off MATLAB:subscripting:noSubscriptsSpecified
end
if ~isempty(dynopt) if ~isempty(dynopt)
% Should be a list of options for the preprocessor in a cell % Should be a list of options for the preprocessor in a cell
...@@ -111,7 +113,11 @@ for i=1:length(varargin) ...@@ -111,7 +113,11 @@ for i=1:length(varargin)
end end
end end
eqlist = eqlist(1:eqnum,:); eqlist = eqlist(1:eqnum,:);
if isoctave && octave_ver_less_than('6')
[~, idx] = unique_stable(eqlist(:,1));
else
[~, idx] = unique(eqlist(:,1), 'stable'); [~, idx] = unique(eqlist(:,1), 'stable');
end
eqlist = eqlist(idx, :); eqlist = eqlist(idx, :);
% Get endogenous variables. % Get endogenous variables.
...@@ -134,7 +140,11 @@ for i=1:length(varargin) ...@@ -134,7 +140,11 @@ for i=1:length(varargin)
fclose(fid); fclose(fid);
end end
elist = elist(1:enum,:); elist = elist(1:enum,:);
if isoctave && octave_ver_less_than('6')
[~, idx] = unique_stable(elist(:,1));
else
[~, idx] = unique(elist(:,1), 'stable'); [~, idx] = unique(elist(:,1), 'stable');
end
elist = elist(idx,:); elist = elist(idx,:);
% Get exogenous variables. % Get exogenous variables.
...@@ -162,7 +172,11 @@ for i=1:length(varargin) ...@@ -162,7 +172,11 @@ for i=1:length(varargin)
fclose(fid); fclose(fid);
end end
xlist = xlist(1:xnum,:); xlist = xlist(1:xnum,:);
if isoctave && octave_ver_less_than('6')
[~, idx] = unique_stable(xlist(:,1));
else
[~, idx] = unique(xlist(:,1), 'stable'); [~, idx] = unique(xlist(:,1), 'stable');
end
xlist = xlist(idx,:); xlist = xlist(idx,:);
% Get parameter values. % Get parameter values.
...@@ -255,8 +269,9 @@ end ...@@ -255,8 +269,9 @@ end
fprintf(fid, 'end;'); fprintf(fid, 'end;');
fclose(fid); fclose(fid);
if ~isoctave && matlab_ver_less_than('9.14')
warning on MATLAB:subscripting:noSubscriptsSpecified warning on MATLAB:subscripting:noSubscriptsSpecified
end
function b = isequationtag(str) function b = isequationtag(str)
b = true; b = true;
......
function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M) function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, T, M)
% Wrapper around the *_static.m file, for use with dynare_solve, % Wrapper around the *_static.m file, for use with dynare_solve,
% when block_mfs option is given to steady. % when block_mfs option is given to steady.
% Copyright (C) 2009-2020 Dynare Team % Copyright (C) 2009-2022 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -21,4 +21,4 @@ function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M) ...@@ -21,4 +21,4 @@ function [r, g1] = block_bytecode_mfs_steadystate(y, b, y_all, exo, params, M)
indx = M.block_structure_stat.block(b).variable; indx = M.block_structure_stat.block(b).variable;
y_all(indx) = y; y_all(indx) = y;
[r, g1] = bytecode( y_all, exo, params, y_all, 1, y_all, 'evaluate', 'static', ['block = ' int2str(b) ]); [r, g1] = bytecode(y_all, exo, params, y_all, 1, y_all, T, 'evaluate', 'static', ['block = ' int2str(b) ]);
...@@ -81,8 +81,8 @@ for draw=1:options_.bvar_replic ...@@ -81,8 +81,8 @@ for draw=1:options_.bvar_replic
% Build the IRFs... % Build the IRFs...
lags_data = zeros(ny,ny*nlags) ; lags_data = zeros(ny,ny*nlags) ;
sampled_irfs(:,:,1,draw) = Sigma_lower_chol ; sampled_irfs(:,:,1,draw) = StructuralMat ;
lags_data(:,1:ny) = Sigma_lower_chol ; lags_data(:,1:ny) = StructuralMat ;
for t=2:options_.irf for t=2:options_.irf
sampled_irfs(:,:,t,draw) = lags_data(:,:)*Phi(1:ny*nlags,:) ; sampled_irfs(:,:,t,draw) = lags_data(:,:)*Phi(1:ny*nlags,:) ;
lags_data(:,ny+1:end) = lags_data(:,1:end-ny) ; lags_data(:,ny+1:end) = lags_data(:,1:end-ny) ;
......
...@@ -66,7 +66,7 @@ Parameters = Parameters'; ...@@ -66,7 +66,7 @@ Parameters = Parameters';
try try
dd = chol(CovarianceMatrix); dd = chol(CovarianceMatrix);
catch catch
error('The covariance matrix has to be a symetric positive definite matrix!') error('The covariance matrix has to be a symmetric positive definite matrix!')
end end
% Set parameters related to the proposal distribution % Set parameters related to the proposal distribution
......
function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M) function estim_params=check_for_calibrated_covariances(estim_params,M_)
% function check_for_calibrated_covariances(xparam1,estim_params,M) % function check_for_calibrated_covariances(estim_params,M)
% find calibrated covariances to consider during estimation % find calibrated covariances to consider during estimation
% Inputs % Inputs
% -xparam1 [vector] parameters to be estimated
% -estim_params [structure] describing parameters to be estimated % -estim_params [structure] describing parameters to be estimated
% -M [structure] describing the model % -M_ [structure] describing the model
% %
% Outputs % Outputs
% -estim_params [structure] describing parameters to be estimated % -estim_params [structure] describing parameters to be estimated
...@@ -12,7 +11,7 @@ function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M) ...@@ -12,7 +11,7 @@ function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M)
% Notes: M is local to this function and not updated when calling % Notes: M is local to this function and not updated when calling
% set_all_parameters % set_all_parameters
% Copyright (C) 2013-2017 Dynare Team % Copyright © 2013-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -28,27 +27,108 @@ function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M) ...@@ -28,27 +27,108 @@ function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M)
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
Sigma_e_calibrated=M.Sigma_e;
H_calibrated=M.H; if isfield(estim_params,'calibrated_covariances')
%check covariance for structural errors estim_params = rmfield(estim_params,'calibrated_covariances'); %remove if already present
covariance_pos=find(tril(Sigma_e_calibrated,-1)); %off-diagonal elements set by covariances before updating correlation matrix to reflect estimated covariances end
covariance_pos_ME=find(tril(H_calibrated,-1)); %off-diagonal elements set by covariances before updating correlation matrix to reflect estimated covariances if isfield(estim_params,'calibrated_covariances_ME')
estim_params = rmfield(estim_params,'calibrated_covariances_ME'); %remove if already present
%locally updated M end
M = set_all_parameters(xparam1,estim_params,M);
[rows_calibrated, columns_calibrated]=ind2sub(size(M_.Sigma_e),find(tril(M_.Sigma_e,-1))); %find linear indices of preset lower triangular covariance entries
correlation_pos=find(tril(M.Correlation_matrix,-1)); %off-diagonal elements set by correlations after accounting for estimation
calibrated_covariance_pos=covariance_pos(~ismember(covariance_pos,correlation_pos)); if estim_params.ncx %delete preset entries actually estimated
if any(calibrated_covariance_pos) for i=1:estim_params.ncx
[rows, columns]=ind2sub(size(M.Sigma_e),calibrated_covariance_pos); %find linear indices of lower triangular covariance entries shock_1 = estim_params.corrx(i,1);
estim_params.calibrated_covariances.position=[calibrated_covariance_pos;sub2ind(size(M.Sigma_e),columns,rows)]; %get linear entries of upper triangular parts shock_2 = estim_params.corrx(i,2);
estim_params.calibrated_covariances.cov_value=Sigma_e_calibrated(estim_params.calibrated_covariances.position); estimated_corr_pos=find(rows_calibrated==shock_1 & columns_calibrated==shock_2);
end if ~isempty(estimated_corr_pos)
rows_calibrated(estimated_corr_pos)=[];
correlation_pos_ME=find(tril(M.Correlation_matrix_ME,-1)); %off-diagonal elements set by correlations after accounting for estimation columns_calibrated(estimated_corr_pos)=[];
calibrated_covariance_pos_ME=covariance_pos_ME(~ismember(covariance_pos_ME,correlation_pos_ME)); end
if any(calibrated_covariance_pos_ME) estimated_corr_pos=find(rows_calibrated==shock_2 & columns_calibrated==shock_1);
[rows, columns]=ind2sub(size(M.H),calibrated_covariance_pos_ME); %find linear indices of lower triangular covariance entries if ~isempty(estimated_corr_pos)
estim_params.calibrated_covariances_ME.position=[calibrated_covariance_pos_ME;sub2ind(size(M.H),columns,rows)]; %get linear entries of upper triangular parts rows_calibrated(estimated_corr_pos)=[];
estim_params.calibrated_covariances_ME.cov_value=H_calibrated(estim_params.calibrated_covariances_ME.position); columns_calibrated(estimated_corr_pos)=[];
end
end
if any(rows_calibrated)
estim_params.calibrated_covariances.position=[sub2ind(size(M_.Sigma_e),rows_calibrated,columns_calibrated);sub2ind(size(M_.Sigma_e),columns_calibrated,rows_calibrated)]; %get linear entries of upper triangular parts
estim_params.calibrated_covariances.cov_value=M_.Sigma_e(estim_params.calibrated_covariances.position);
end
end
[rows_calibrated, columns_calibrated]=ind2sub(size(M_.H),find(tril(M_.H,-1))); %find linear indices of preset lower triangular covariance entries
if estim_params.ncn %delete preset entries actually estimated
for i=1:estim_params.ncn
shock_1 = estim_params.corrn(i,1);
shock_2 = estim_params.corrn(i,2);
estimated_corr_pos=find(rows_calibrated==shock_1 & columns_calibrated==shock_2);
if ~isempty(estimated_corr_pos)
rows_calibrated(estimated_corr_pos)=[];
columns_calibrated(estimated_corr_pos)=[];
end
estimated_corr_pos=find(rows_calibrated==shock_2 & columns_calibrated==shock_1);
if ~isempty(estimated_corr_pos)
rows_calibrated(estimated_corr_pos)=[];
columns_calibrated(estimated_corr_pos)=[];
end
end
end
if any(rows_calibrated)
estim_params.calibrated_covariances_ME.position=[sub2ind(size(M_.H),rows_calibrated,columns_calibrated);sub2ind(size(M_.H),columns_calibrated,rows_calibrated)]; %get linear entries of upper triangular parts
estim_params.calibrated_covariances_ME.cov_value=M_.H(estim_params.calibrated_covariances_ME.position);
end
return % --*-- Unit tests --*--
%@test:1
M_.Sigma_e=[1 0; 0 1];
M_.H=[1 0; 0 1];
M_.Correlation_matrix= [1 -0.5; -0.5 1];
M_.Correlation_matrix_ME=[1 -0.5; -0.5 1];
estim_params.ncx=1;
estim_params.ncn=1;
estim_params.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params.corrn=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params=check_for_calibrated_covariances(estim_params,M_);
if isfield(estim_params,'calibrated_covariances_ME') || isfield(estim_params,'calibrated_covariances')
t(1)=false;
else
t(1)=true;
end
M_.Sigma_e=[1 -0.1; -0.1 1];
M_.H=[1 -0.1; -0.1 1];
M_.Correlation_matrix= [1 -0.5; -0.5 1];
M_.Correlation_matrix_ME=[1 0; 0 1];
estim_params.ncx=1;
estim_params.ncn=0;
estim_params.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params.corrn=[];
estim_params=check_for_calibrated_covariances(estim_params,M_);
t(2)=isequal(estim_params.calibrated_covariances_ME.position,[2;3]);
t(3)=isequal(estim_params.calibrated_covariances_ME.cov_value,[-0.1;-0.1]);
M_.Sigma_e=[1 -0.1; -0.1 1];
M_.H=[1 -0.1; -0.1 1];
M_.Correlation_matrix= [1 -0.5; -0.5 1];
M_.Correlation_matrix_ME=[1 0; 0 1];
estim_params.ncx=1;
estim_params.ncn=1;
estim_params.corrx=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params.corrn=[2 1 NaN -1 1 3 0 0.2000 NaN NaN NaN];
estim_params=check_for_calibrated_covariances(estim_params,M_);
if isfield(estim_params,'calibrated_covariances_ME') || isfield(estim_params,'calibrated_covariances')
t(4)=false;
else
t(4)=true;
end end
T = all(t);
%@eof:1